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

[Бахурин Сергей]

Промоделировал все методы из статьи. Имеется небольшая ошибка на выходе. Ищу причину. Визуально все должно быть ок. По итогам отпишусь.

[Алексей Климов]

Есть еще вот эта статья:
An Enhanced Interpolated-Modulated Sliding DFT for High Reporting Rate PMUs
Paolo Romano, Mario Paolone

Моя моделька почти повторяет рисунок авторов:

И они же исследуют стабильность, которая явно лучше моей:

Хотя, они молчат о типе данных...

[Бахурин Сергей]

Проверил три варианта:
1) Классическое скользящее ДПФ

1.png (5.51 КБ) 7148 просмотров

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

2.png (6.69 КБ) 7148 просмотров

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

3.png (5.26 КБ) 7148 просмотров

Код симуляции (делал в 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)');

В результате все три варианта возвращают малую ошибку:

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

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

cic.png (4.77 КБ) 7148 просмотров

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

[Алексей Климов]

Спасибо за такое внимание к моей проблеме.

Возможно, вы можете подсказать где дальше копать информацию о целочисленной реализации (как посчитать количество требуемых разрядов и прочие примитивные вопросы)?

[Бахурин Сергей]

Проблема в том, что операции с плавающей точкой неизбежно вносят ошибки округления, даже суммирование вносит небольшую ошибку на уровне представления числа с плавающей точкой. Поэтому интегратор копит эту ошибку и через некоторое время она неизбежно станет большой. Но можно исключить накопление ошибок если свести расчет к целочисленной арифметике.

Добавьте к скрипту выше ещё один вариант

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

% Варинт 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)');

Можно видеть что ошибка всегда есть, но ее уровень не растет, потому что операции без округления. Уровень ошибки зависит теперь не от ошибок округления вычислителя, а от ошибок перехода к целочисленной арифметике (парметра G )

[Алексей Климов]

cic.PNG (7.26 КБ) 7141 просмотр

Что-то такое ведь? Accumulator сумматора тоже в int32. Но неустойчивость проявляется опять

[Бахурин Сергей]

Вот за это я не люблю всякие симулинки. Как преобразование делается к инт32? Если через простой floor то неустойчивость конечно проявится. Поставьте перед интегратором гребенчатой звено с вычитанием через задержку N. По идее постоянная составляющая исчезнет и должно всё быть устойчиво.

[Алексей Климов]

Я попробовал все, изначально выбирал zero (несмещенный):

mode.png (4.64 КБ) 7138 просмотров

[Бахурин Сергей]

Может до округления была постоянная составляющая на входе?

[Алексей Климов]

В итоге моделька заработала, видимо где-то накосячил, но не заметил.
Преобразование к инт32 сделал c параметром zero после умножения на 2^20 (так получилось даже не сильно увеличить затраты в плис).
Благодарю.