【问题标题】:Comparing Naive Inverse Filter to Wiener Filter for Deconvolution in Matlab在 Matlab 中比较朴素逆滤波器和维纳滤波器以进行反卷积
【发布时间】:2014-12-06 10:22:46
【问题描述】:

我目前正在尝试将一个简单的逆滤波器与维纳滤波器进行比较,以使用 matlab 进行反卷积。我的起始信号是exp(-t^2),这将与一个非零的矩形进行卷积,时间为-.5 到 0.5。我正在引入幅度在 -.5 到 .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

对于第一种方法,我只是将 c 的 FT 除以 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');

【问题讨论】:

  • 去除噪声并计算简单的逆滤波器会显示原始矩形(前提是我已经更改了上面的代码以反映)。我希望简单逆滤波器的输出能够给出较大的值,因为我主要是除以较小的值。这或多或少与我实际得到的输出相匹配。我认为我现在的主要问题是计算维纳滤波器。我已经更新了我的代码以反映我现在的想法,但我对此非常不确定,它仍然不会产生任何像原始矩形一样的东西。
  • 我还尝试通过直接计算我认为是 I 和 R 的两侧功率谱密度来计算维纳滤波器。我更新了上面的代码以反映这一点。我现在得到类似 sinc 的东西。所以这可能会更好,但仍然没有关闭。
  • 我也尝试过两次关于维纳滤波的 ifft,我知道这很愚蠢,但它给出了正确幅度的矩形(因为第一个 ifft 是 sinc),但是宽度是错误的...我已经更新了上面的代码以显示此行。

标签: matlab signal-processing fft convolution ifft


【解决方案1】:

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

【讨论】:

  • 我试图做同样的事情,但您的代码现在除以零,导致所有 NaN 用于维纳反卷积。
  • 是的,您必须检查 NaN 值。
  • 您好,我尝试运行您的代码,但它说 I 和 R 未定义。 I 和 R 的值是多少?
  • @MoneyBall 我上次看这个已经有一段时间了,但我相信 I = FxS 其中 F = fft(f) 和 S = fft(s)(f是原始信号,s是矩形)
猜你喜欢
  • 1970-01-01
  • 1970-01-01
  • 2017-11-22
  • 2021-01-03
  • 2017-11-04
  • 2019-11-24
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
相关资源
最近更新 更多