Доброго времени суток. Могли бы сказать что за 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
который можно сделать гарантированно устойчивым, если выполнять в целочисленной арифметике. В этом случае необходимо иметь на входе целочисленный сигнал после умножителя, и выполнять все сложения без округлений. Тогда система будет устойчивой сколь угодно долго.
Скользящие ДПФ
Re: Скользящие ДПФ
Re: Скользящие ДПФ
Пройдите вверх по ветке и там будет архив doc c этим файлом. Но там просто запись тестового сигнал, можете сгенерировать свой.
- Бахурин Сергей
- Администратор
- Сообщения: 1116
- Зарегистрирован: 05 окт 2010, 19:55
- Контактная информация:
-
- Сообщения: 1
- Зарегистрирован: 26 мар 2021, 12:51
Re: Скользящие ДПФ
Сергей, можете уточнить, пожалуйста. Вы указали, что требуется коррекция по фазе на выходе. Что Вы имели в виду? И почему возникает такая необходимость?Бахурин Сергей писал(а): ↑17 дек 2018, 12:46Вопрос такой. Насколько важно корректность фазы? Или амплитуды будет достаточно?
Ещё вопрос в схеме задержка , чему равно d?
Попробовал ваши данные. Если все в формате double, то все устойчиво. Неустойчивость возможна если сильно округлить поворотный коэффициент (до 4 знаков после запятой в большую сторону). Но замечу что задержка d =M+1, а не M (часто сам ошибаюсь).
Даже в этом случае можно пытаться стабилизировать притягивая полюс на ед. окружность.
Можно попробовать совсем исключить потенциальную неустойчивость путем изменения структуры на ту что приведена в статье, но возникает необходимость коррекции фазы на выходе. Если вам достаточно только амплитуды, то это будет проще.