fiona_chen 发表于 2007-5-30 16:03

理想低通滤波器的实现

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

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

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

鞠躬拜谢了!

recbio 发表于 2007-5-31 00:23

这个问题源于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

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

fiona_chen 发表于 2007-5-31 03:51

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

freecat 发表于 2007-6-22 22:56

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

freecat 发表于 2007-6-22 23:03

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

recbio 发表于 2007-6-23 22:02

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

首先我们说一下这个中间的桥梁——拉氏变换——就是用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]
查看完整版本: 理想低通滤波器的实现