Дискретизация и восстановление сигнала в octave. Теорема Котельникова
Продолжение серии заметок по работе в octave. В этот раз это цифровая обработка сигналов и теорема Котельникова, которая гласит, что если максимальная частота в сигнале больше половины частоты его дискретизации, то сигнал не может быть восстановлен без искажений. В следующем коде это утверждение проверяется, применяя различную частоту дискретизации.
set (0, 'defaulttextfontname', 'Terminus,16');
QUALITY = 16;
N = 2048; % Максимальное значение времени t
freq = 257/(N-1); % Частота волны
dfreq = 2 * freq; % Частота дискретизации (замеров)
disp('Частота волны freq = '), disp(freq)
t = 0:1/(20*dfreq):N;
% Уравнение волны
F = N*log(N+1)/(N^2-1)*sin((2*pi*100/N)*t+N/10) + \
(N+3)*log(N/3)/(N^2-1)*sin((2*pi*157/N)*t+N/3) + \
(N+4)*log(N/4)/(N^2-1)*cos((2*pi*257/(N-1))*t+N/4);
% Уравнение волны 0.7
t1 = 0:1/(0.7*dfreq):N;
F1 = N*log(N+1)/(N^2-1)*sin((2*pi*100/N)*t1+N/10) + \
(N+3)*log(N/3)/(N^2-1)*sin((2*pi*157/N)*t1+N/3) + \
(N+4)*log(N/4)/(N^2-1)*cos((2*pi*257/(N-1))*t1+N/4);
t2 = 0:1/(1*dfreq):N;
% Уравнение волны 1
F2 = N*log(N+1)/(N^2-1)*sin((2*pi*100/N)*t2+N/10) + \
(N+3)*log(N/3)/(N^2-1)*sin((2*pi*157/N)*t2+N/3) + \
(N+4)*log(N/4)/(N^2-1)*cos((2*pi*257/(N-1))*t2+N/4);
t3 = 0:1/(1.5*dfreq):N;
% Уравнение волны 1.5
F3 = N*log(N+1)/(N^2-1)*sin((2*pi*100/N)*t3+N/10) + \
(N+3)*log(N/3)/(N^2-1)*sin((2*pi*157/N)*t3+N/3) + \
(N+4)*log(N/4)/(N^2-1)*cos((2*pi*257/(N-1))*t3+N/4);
% Уравнение волны: A * sin (2 \pi t\ omega)
% Имеем гармонику с наибольшей частотой 2*pi*257/(N-1)*t: частота = 257/(n-1)
% По т. Котельникова частота дискретизации должна быть в 2 раза больше частоты самой высокочастотной гармоники, поэтому она равна 2 \omega
figure(1) % Исходный сигнал
subplot(3,1,1)
plot(t1,F1)
axis([0. 50])
subplot(3,1,2)
plot(t2,F2)
axis([0. 50])
subplot(3,1,3)
plot(t3,F3)
axis([0. 50])
% Интерполяция 0.7
figure(2)
subplot(3,1,1);
plot(t1, F1)
axis([0. 50])
subplot(3,1,2);
ti1 = 0:1/(0.7*dfreq*1.5):N;
FI1 = interp1(t1,F1,ti1,'spline');
plot(ti1, FI1)
axis([0. 50])
subplot(3,1,3);
plot(t,F)
axis([0. 50])
% Интерполяция 1
figure(3)
subplot(3,1,1);
plot(t2, F2)
axis([0. 50])
subplot(3,1,2);
ti2 = 0:1/(1*dfreq*1.5):N;
FI2 = interp1(t2,F2,ti2,'spline');
plot(ti2, FI2)
axis([0. 50])
subplot(3,1,3);
plot(t,F)
axis([0. 50])
% Интерполяция 1.5
figure(4)
subplot(3,1,1);
plot(t3, F3)
axis([0. 50])
subplot(3,1,2);
ti3 = 0:1/(1.5*dfreq*1.5):N;
FI3 = interp1(t3,F3,ti3,'spline');
plot(ti3, FI3)
axis([0. 50])
subplot(3,1,3);
plot(t,F)
axis([0. 50])
Ответить