Re: Скользящие ДПФ
Добавлено: 18 дек 2018, 12:14
Промоделировал все методы из статьи. Имеется небольшая ошибка на выходе. Ищу причину. Визуально все должно быть ок. По итогам отпишусь.
Теория и практика цифровой обработки сигналов
http://ru.dsplib.org/forum/
Код: Выделить всё
clear all; clc;
%graphics_toolkit("gnuplot");
load data.mat
M = 1152
N = length(x);
k = 3;
P = 20000;
% Скользящее дпф в лоб
for m = 1:P-M+10
tmp = fft(x(m:m+M-1));
Z(m+M-1) = tmp(k+1);
end
w = exp(-2i*pi*k / M);
b = zeros(M+1, 1);
b(1) = 1;
b(M+1) = -1;
a = [w, -1];
Y = filter(b,a,x);
figure(1); subplot(411); plot(1:P, real(Z(1:P)), 1:P, real(Y(1:P)));
title('Классический вариант');
ylabel('Re(Z), Re(Y)');
figure(1); subplot(412); plot(M+1:P, real(Z(M+1:P).' - Y(M+1:P)));
ylabel('Re(Z-Y)');
figure(1); subplot(413); plot(1:P, imag(Z(1:P)), 1:P, imag(Y(1:P)));
ylabel('Im(Z), Im(Y)');
figure(1); subplot(414); plot(M+1:P, imag(Z(M+1:P).' - Y(M+1:P)));
ylabel('Im(Z-Y)');
% Варинт 1
t = (0:N-1)';
y = x.*exp(-2i * pi * t * k / M);
b = zeros(M+1, 1)';
b(1) = 1;
b(M+1) = -1;
a = [1; -1];
v = filter(b,a,y);
Y = v.*exp(2i*pi* k* (t+1)/M);
figure(2); subplot(411); plot(1:P, real(Z(1:P)), 1:P, real(Y(1:P)));
title('Вариант 1');
ylabel('Re(Z), Re(Y)');
figure(2); subplot(412); plot(M+1:P, real(Z(M+1:P).' - Y(M+1:P)));
ylabel('Re(Z-Y)');
figure(2); subplot(413); plot(1:P, imag(Z(1:P)), 1:P, imag(Y(1:P)));
ylabel('Im(Z), Im(Y)');
figure(2); subplot(414); plot(M+1:P, imag(Z(M+1:P).' - Y(M+1:P)));
ylabel('Im(Z-Y)');
% Варинт 2
t = (0:N-1)';
y = x.*exp(-2i * pi * t * k / M);
b = zeros(M+1, 1)';
b(1) = 1;
b(M+1) = -1;
a = [1; -1];
y = filter(b,1,x);
y = y.*exp(-2i*pi* k* t/M);
v = filter(1,a,y);
Y = v.*exp(2i*pi* k* (t+1)/M);
figure(3); subplot(411); plot(1:P, real(Z(1:P)), 1:P, real(Y(1:P)));
title('Вариант 2');
ylabel('Re(Z), Re(Y)');
figure(3); subplot(412); plot(M+1:P, real(Z(M+1:P).' - Y(M+1:P)));
ylabel('Re(Z-Y)');
figure(3); subplot(413); plot(1:P, imag(Z(1:P)), 1:P, imag(Y(1:P)));
ylabel('Im(Z), Im(Y)');
figure(3); subplot(414); plot(M+1:P, imag(Z(M+1:P).' - Y(M+1:P)));
ylabel('Im(Z-Y)');
Код: Выделить всё
% Варинт 4
t = (0:N-1)';
G=1000000;
y =round(G* x.*exp(-2i * pi * t * k / M));
v = x*0;
Y = x*0;
v(1)=y(1);
Y(1) = v(1);
for m=2:N
if m <M+1
v(m) = y(m);
else
v(m)=y(m) - y(m-M);
end
Y(m) = Y(m-1) + v(m);
end
Y = v.*exp(2i*pi* k* (t+1)/M)/G;
figure(4); subplot(411); plot(1:P, real(Z(1:P)), 1:P, real(Y(1:P)));
title('Вариант 1');
ylabel('Re(Z), Re(Y)');
figure(4); subplot(412); plot(M+1:P, real(Z(M+1:P).' - Y(M+1:P)));
ylabel('Re(Z-Y)');
figure(4); subplot(413); plot(1:P, imag(Z(1:P)), 1:P, imag(Y(1:P)));
ylabel('Im(Z), Im(Y)');
figure(2); subplot(414); plot(M+1:P, imag(Z(M+1:P).' - Y(M+1:P)));
ylabel('Im(Z-Y)');