Re: Скользящие ДПФ
Добавлено: 09 дек 2020, 04:13
Доброго времени суток. Могли бы сказать что за load data.mat в начале кода. Попытался запустить на Octave и не выходит.Бахурин Сергей писал(а): ↑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
который можно сделать гарантированно устойчивым, если выполнять в целочисленной арифметике. В этом случае необходимо иметь на входе целочисленный сигнал после умножителя, и выполнять все сложения без округлений. Тогда система будет устойчивой сколь угодно долго.