Сравнение наивного обратного фильтра с фильтром Винера для деконволюции в Matlab

В настоящее время я пытаюсь сравнить простой обратный фильтр с фильтром Винера для деконволюции с использованием Matlab. Мой начальный сигнал - exp(-t^2), и его нужно свернуть с прямоугольником, который не равен нулю для времени от -.5 до .5. Я ввожу шум с амплитудой в диапазоне от -0,5 до 0,5.

Определение моей временной области для сопоставления частотной области:

f = exp(-t^2) => F

s = rect => R

c = f*s => C

r = noise (see above) => R

with noise c becomes: c = f*s + n => C = FxS + N

Для первого подхода я просто беру FT c и делю его на FT f, а затем выполняю обратный FT. Это составляет s = (approx.) ifft((FxS + N)/F).

Для второго подхода я беру фильтр Винера W и умножаю его на C/R, а затем выполняю обратный FT. Это составляет S = (approx.) ifft(CxW/R).

Винеровский фильтр W = mag_squared(FxS)/(mag_squared(FxS) + mag_squared(N)).

Я использовал «*» для обозначения свертки и «x» для обозначения умножения.

Я пытаюсь сравнить две деконволюции прямоугольника во временном интервале от -3 до 3. Прямо сейчас мои полученные графики деконволюционного прямоугольника совсем не похожи на оригинал.
Может ли кто-нибудь указать мне правильное направление относительно того, что Я делаю неправильно? Я пробовал использовать ifftshift и разные масштабы во многих разных порядках, но, похоже, ничего не работает.

Спасибо

Мой код Matlab ниже:

%%using simple inverse filter
dt = 1/1000;
t = linspace(-3,3,1/dt); %time
s = zeros(1,length(t)); 
s(t>=-0.5 & t<=0.5) = 1; %rect
f = exp(-(t.^2)); %function
r = -.5 + rand(1,length(t)); %noise

S = fft(s);
F = fft(f);
R = fft(r);
C = F.*S + R;
S_temp = C./F;
s_recovered_1 = real(ifft(S_temp));  %correct?...works for signal without R (noise)

figure();
plot(t,s + r);
title('rect plus noise');

figure();
hold on;
plot(t,s,'r');
plot(t,f,'b');
legend('rect input','function');
title('inpute rect and exponential functions');
hold off;

figure();
plot(t,s_recovered_1,'black');
legend('recovered rect');
title('recovered rect using naive filter');


%% using wiener filter
N = length(s);
I_mag = abs(I).^2;
R_mag = abs(R).^2;
W = I_mag./(I_mag + R_mag);
S_temp = (C.*W)./F;
s_recovered_2 = abs(ifft(S_temp));  

figure();
freq = -fs/2:fs/N:fs/2 - fs/N;
hold on;
plot(freq,10*log10(I_mag),'r');
plot(freq,10*log10(R_mag),'b');
grid on
legend('I_mag','R_mag');
title('Periodogram Using FFT')
xlabel('Frequency (Hz)')
ylabel('Power/Frequency (dB/Hz)')

figure();
plot(t,s_recovered_2);
legend('recovered rect');
title('recovered rect using wiener filter');

person John Cast    schedule 11.10.2014    source источник
comment
Удаление шума и вычисление простого обратного фильтра показывают исходный прямоугольник (при условии, что я делаю 's_recovered_1 = real(ifft(S_temp));', что я изменил в приведенном выше коде, чтобы отразить его). Я ожидаю, что на выходе простого обратного фильтра будут большие значения, так как я в основном делю на маленькие значения. Это более или менее соответствует результату, который я действительно получаю. Я думаю, что моя главная проблема сейчас заключается в вычислении фильтра Винера. Я обновил свой код, чтобы отразить то, что я думаю сейчас, но я очень не уверен в этом, и он по-прежнему не создает ничего похожего на исходный прямоугольник.   -  person John Cast    schedule 12.10.2014
comment
Я также попытался вычислить фильтр Винера, напрямую вычислив то, что, как я думаю, является двусторонней спектральной плотностью мощности I и R. Я обновил приведенный выше код, чтобы отразить это. Я получаю что-то вроде sinc сейчас. Так что это может быть лучше, но это все еще не близко.   -  person John Cast    schedule 12.10.2014
comment
и я также пытался дважды взять ifft в отношении винеровской фильтрации, что, я знаю, глупо, но дает прямоугольник правильной амплитуды (поскольку первый ifft - это sinc), но ширина неверна... я обновили мой код выше, чтобы показать строку для этого.   -  person John Cast    schedule 12.10.2014


Ответы (1)


Вот и получается, что я делил не на тот знаменатель при расчете фильтра Винера. Теперь я также рассчитываю |...|^2 (спектральную плотность мощности) каждого члена в фильтре Винера, используя простой способ abs(...)^2. Приведенный выше код отражает эти изменения. Надеюсь, это будет полезно для таких нубов, как я :)

person John Cast    schedule 13.10.2014
comment
Я пытаюсь сделать то же самое, но ваш код теперь делится на ноль, что приводит ко всем NaN для деконволюции Винера. - person Leo; 04.11.2014
comment
Да, вы должны проверить значения NaN. - person John Cast; 05.11.2014
comment
Привет, я попытался запустить ваш код, но он говорит, что I и R не определены. Каковы значения I и R? - person MoneyBall; 23.09.2017
comment
@MoneyBall Прошло много времени с тех пор, как я в последний раз смотрел на это, но я считаю, что I = FxS, где F = fft(f) и S = fft(s) (f — исходный сигнал, s — прямоугольник) - person John Cast; 25.09.2017