理想低通滤波器的实现

fiona_chen

研究的头都大了。哪位高手来帮个忙吧。
我现在要实现一个理想低通滤波器，简单说就是频域上的一个门函数，在|Ω|<Ωg的时候是1，sonst（im Bereich Ωg<|Ω|<pi）的时候是0

现在希望在matlab里面实现，并且画出图像，关键是用fft（或者用ifft）求出时域。
或者如果输入时域的sinc函数，然后用fft变到频域也行。
我之前时尝试从时域变到频域，题目说可以用2000点的fft，但是变出来死活不是一个门函数。。。唉～～

帮帮忙吧。看到Skript有类似的，但是就给了图，没说怎么画的。

鞠躬拜谢了！

recbio

这个问题源于matlab一维fft的算法。由于用了蝶形算法，但是没有做相应的对称修正，所以matlab的fft和实际得到结果并不相同。
最简单的检验：
x=0:1999;
y1=sin(x./10*2*pi);    (t=2000/10 = 100 假设的高频)
y2=sin(x/100*2*pi);   (t=2000/100= 10 假设的低频)
y=y1+y2;
fy=fft(y);
plot(abs(fy));
你看到了什么？是四个峰！并且，我还可以告诉你，中间的两个是高频的峰，外面的两个是低频的峰。且峰的位置还比较怪！
你怎么用单边的函数来滤波？

一个简单的方法是，只用一半。不过，会引入失真，尤其是低通；
一个比较保守的方法，就是用fftshft来修正；
如：
figure;plot(fftshift(fy));
你会看到，中间的两个峰是低频峰，外面的两个峰是高频峰。只需要从中间开始cut，你就可以除去高频了；
filter=[repmat(0,1,900) repmat(1,1,200) repmat(0,1,900)] ;
中间有200个1 转折在  (1000-900)/2 = 50。既<50都留下；
yf=real(ifft(fftshift((fftshift(fy).*filter))));
再看看yf；
figure;plot(fftshift(yf));
怎么样，是不是和y2一样？看看误差是多少：
sqrt(sum((yf-y2).^2)/2000)

ans =

  1.3664e-014

这个值基本上可以认可了。

fiona_chen

谢谢ls的。程序终于有点样子了。

freecat

其实不用fft,ifft计算也可以
直接从频率域z变换到时间域即可得解

freecat

接下来就可以用各种窗函数造各种低通了~

recbio

楼上的建议非常好。其实，如果知道这些变换的关系，对于信号处理很有用。

首先我们说一下这个中间的桥梁——拉氏变换——就是用exp(-sT)来卷积原函数。
其中s ＝σ＋jΩ ，T是周期。j^2=-1；σ 和 Ω 是实数。

于是，我们看到，假如 σ=0 这个变换就是 付氏变换 了。

那么  拉氏变换 和 z变换之间的关系呢？
其实，对于函数的z变换，假如 z= exp(-sT)，那么，实数在z变换平面的对应关系，等于实数在拉氏变换平面的对应关系。
同理，我们很容易发现，在这个时候，付氏变换就是z平面的一个圆了，
因为拉氏变换中 s ＝σ＋jΩ 而对应的  z变换    z ＝r*exp(jω)
如果是付氏变换 σ=0 ，  令 z= exp(-sT)  解出 r＝1

可以说， 付氏变换 是 z变换平面 在 单位圆 里面的 拉氏变换。

当然，解一般低通滤波用 z变换 就象用牛刀来杀鸡，一定是“ 游刃有余”了。