Скользящие ДПФ

Олег
Сообщения: 1
Зарегистрирован: 09 дек 2020, 04:08

Re: Скользящие ДПФ

Сообщение Олег »

Бахурин Сергей писал(а):
18 дек 2018, 20:13
Проверил три варианта:
1) Классическое скользящее ДПФ
1.png

2) Вариант 1 с выносом умножителя наружу
2.png


3) Варинат два с размещением умножителя между дифференциальной частью и интегратором
3.png


Код симуляции (делал в octave, в матлабе должно заработать

Код: Выделить всё

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)');
В результате все три варианта возвращают малую ошибку:
plot.png

Причем вариант 1 и 2 дают на порядок меньшую ошибку классического варианта. Но во всех трех случаях наблюдается рост ошибки, вызванный накоплением ошибок округления в рекурсивной части фильтра. Я не смог получить результата, когда ошибка варианта 1 и 2 становится сколь угодно весомой, но думаю, что если фильтр не будет сбрасываться довольно долго, то неустойчивость проявит себя.

Однако если взглянуть на вариант 1, то можно видеть, что в нем есть CIC фильтр,

cic.png

который можно сделать гарантированно устойчивым, если выполнять в целочисленной арифметике. В этом случае необходимо иметь на входе целочисленный сигнал после умножителя, и выполнять все сложения без округлений. Тогда система будет устойчивой сколь угодно долго.
Доброго времени суток. Могли бы сказать что за load data.mat в начале кода. Попытался запустить на Octave и не выходит.

kaa
Сообщения: 33
Зарегистрирован: 17 мар 2019, 20:03

Re: Скользящие ДПФ

Сообщение kaa »

Пройдите вверх по ветке и там будет архив doc c этим файлом. Но там просто запись тестового сигнал, можете сгенерировать свой.

Аватара пользователя
Бахурин Сергей
Администратор
Сообщения: 1004
Зарегистрирован: 05 окт 2010, 19:55
Контактная информация:

Re: Скользящие ДПФ

Сообщение Бахурин Сергей »

kaa писал(а):
09 дек 2020, 15:40
Пройдите вверх по ветке и там будет архив doc c этим файлом. Но там просто запись тестового сигнал, можете сгенерировать свой.
Если быть точными, сообщение от Алексей Климов 17 дек 2018, 10:23

Ответить