理想低通滤波器的实现
研究的头都大了。哪位高手来帮个忙吧。我现在要实现一个理想低通滤波器,简单说就是频域上的一个门函数,在|Ω|<Ωg的时候是1,sonst(im Bereich Ωg<|Ω|<pi)的时候是0
现在希望在matlab里面实现,并且画出图像,关键是用fft(或者用ifft)求出时域。
或者如果输入时域的sinc函数,然后用fft变到频域也行。
我之前时尝试从时域变到频域,题目说可以用2000点的fft,但是变出来死活不是一个门函数。。。唉~~
帮帮忙吧。看到Skript有类似的,但是就给了图,没说怎么画的。
鞠躬拜谢了! 这个问题源于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= ;
中间有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
这个值基本上可以认可了。 谢谢ls的。程序终于有点样子了。 其实不用fft,ifft计算也可以
直接从频率域z变换到时间域即可得解 接下来就可以用各种窗函数造各种低通了~ 楼上的建议非常好。其实,如果知道这些变换的关系,对于信号处理很有用。
首先我们说一下这个中间的桥梁——拉氏变换——就是用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变换 就象用牛刀来杀鸡,一定是“ 游刃有余”了。
页:
[1]