在Matlab中将朴素逆滤波器与维纳滤波器进行反卷积的比较 [英] Comparing Naive Inverse Filter to Wiener Filter for Deconvolution in Matlab
问题描述
我目前正在尝试将简单的逆滤波器与维纳滤波器进行比较，以使用Matlab进行反卷积.我的起始信号是 exp(-t ^ 2)，并且将其与-0.5至.5时间不为零的rect进行卷积.我正在引入幅度在-0.5至.5之间的噪声.
I am currently trying to compare a simple inverse filter to the wiener filter for deconvolution using matlab. My starting signal is exp(-t^2) and this is to be convolved with a rect that is nonzero for times -.5 to .5. I am introducing noise with amplitude in the range -.5 to .5.
定义我的时域到频域映射:
Defining my time domain to frequency domain mapping:
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
对于第一种方法，我只是采用 c 的FT并将其除以 f 的FT，然后进行逆FT.这等于 s =(大约)ifft((FxS + N)/F).
For the first approach I am simply taking the FT of c and dividing it by the FT of f and then doing an inverse FT. This is amounts to s = (approx.) ifft((FxS + N)/F).
对于第二种方法，我采用维纳滤波器，即 W ，并将其乘以 C/R ，然后进行逆FT.这等于 S =(大约)ifft(CxW/R).
For the second approach I am taking the wiener filter, W, and multiplying it against C/R and then doing an inverse FT. This amounts to S = (approx.) ifft(CxW/R).
维纳滤波器为 W = mag_squared(FxS)/(mag_squared(FxS)+ mag_squared(N)).
我用'*'表示卷积，用'x'表示乘法.
I have used '*' to mean convolution and 'x' to mean multiplication.
我正在尝试比较在时间间隔-3到3之间rect的两个反卷积.现在，我得到的反卷积rect的图形看起来与原始图形完全不同.
有人能指出我对错的正确方向吗?我尝试以许多不同的顺序使用ifftshift和不同的缩放比例，但似乎没有任何效果.
I am trying to compare the two deconvolutions of the rect over the time interval -3 to 3.
Right now, my resulting graphs of the deconvolved rect look nothing like the original.
Could someone point me in the right direction as to what I'm doing to wrong? I have tried using ifftshift and different scalings in many different orderings but nothing seems to work.
谢谢
我的matlab代码如下:
My matlab code is below:
%%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');
推荐答案
因此，事实证明，在计算维纳滤波器时，我用错误的分母进行了除法.现在，我还使用简单的abs(...)^ 2方法来计算Wiener滤波器中每个项的| ... | ^ 2(功率谱密度).上面的代码反映了这些更改.希望这对像我这样的菜鸟有帮助:)
So it turns out that I was dividing by the wrong denominator when calculating the Wiener Filter. I also now calculate the |...|^2 (power spectral density) of each term in the Wiener Filter using the straightforward abs(...)^2 way. The code above reflects these changes. Hope this is helpful to any noobs like myself :)
这篇关于在Matlab中将朴素逆滤波器与维纳滤波器进行反卷积的比较的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
