2014-10-11 92 views
3

我目前正在尝试使用matlab将简单逆滤波器与维纳滤波器进行反卷积比较。我的起始信号是exp(-t^2),这是用一个不等于-5的时间长度来卷积的。我正在引入振幅范围为-5到.5的噪音。将Naive逆滤波器与Wiener滤波器进行比较,以便在Matlab中解卷积

定义我的时域到频域映射:

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并通过f的FT划分它,然后做逆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'); 
+0

去除噪音和计算简单逆滤波器揭示了原始矩形(假设我做's_recovered_1 =真实(ifft(S_temp));'我已经改变了上述代码来反映)。我期望简单逆滤波器的输出能够给出很大的值,因为我大部分被小的值所分割。这或多或少地匹配我实际得到的输出。我认为我现在的主要问题是计算维纳滤波器。我已更新我的代码以反映我现在的想法,但我对此非常不确定,它仍然不会产生任何类似原始矩形的东西。 – 2014-10-11 22:29:56

+1

我也尝试通过直接计算我认为是I和R的双侧功率谱密度来计算维纳滤波器。我更新了上面的代码以反映这一点。我现在得到了一个像sinc的东西。所以这可能会更好,但它还没有结束。 – 2014-10-11 23:42:49

+0

我也尝试过两次关于维纳滤波,这是愚蠢的,我知道,但它给出了正确的幅度矩(因为第一个ifft是一个sinc),但宽度是错误的... 我已经更新了上面的代码以显示此行。 – 2014-10-11 23:48:50

回答

2

所以,事实证明,我被错误的分母计算维纳滤波器时分裂。我现在也使用简单的abs(...)^ 2方法计算Wiener滤波器中每个项的| ... |^2(功率谱密度)。上面的代码反映了这些变化。 希望这有助于像我这样的小菜鸟:)

+0

即时尝试做同样的事情,但你的代码现在分为零,导致所有NaN的维纳解卷积。 – Leo 2014-11-04 16:55:49

+1

是的,你必须检查NaN值。 – 2014-11-05 12:13:40

+0

嗨,我试着运行你的代码,但它说我和R没有定义。 I和R的价值是什么? – MoneyBall 2017-09-23 08:55:50