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

Олег
Сообщения: 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
Сообщения: 36
Зарегистрирован: 17 мар 2019, 20:03

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

Сообщение kaa »

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

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

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

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

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

Nina Chechulina
Сообщения: 1
Зарегистрирован: 26 мар 2021, 12:51

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

Сообщение Nina Chechulina »

Бахурин Сергей писал(а):
17 дек 2018, 12:46
Вопрос такой. Насколько важно корректность фазы? Или амплитуды будет достаточно?

Ещё вопрос в схеме задержка , чему равно d?

Попробовал ваши данные. Если все в формате double, то все устойчиво. Неустойчивость возможна если сильно округлить поворотный коэффициент (до 4 знаков после запятой в большую сторону). Но замечу что задержка d =M+1, а не M (часто сам ошибаюсь).

Даже в этом случае можно пытаться стабилизировать притягивая полюс на ед. окружность.

Можно попробовать совсем исключить потенциальную неустойчивость путем изменения структуры на ту что приведена в статье, но возникает необходимость коррекции фазы на выходе. Если вам достаточно только амплитуды, то это будет проще.
Сергей, можете уточнить, пожалуйста. Вы указали, что требуется коррекция по фазе на выходе. Что Вы имели в виду? И почему возникает такая необходимость?

Ответить