傅里叶变化与频率那些事儿
学数字信号处理的时候都会接触到的三个东西就是数字频率,模拟频率,采样频率,这三个东西其实挺好懂得,但是一旦涉及到编程那些东西好像很容易出错,而且一般都是比较难查的。然后一旦到了数字滤波器那些东西,又来什么归一化频率的,就更乱了,加上数字福利叶变换后频率就变成-π到π了,然后又有人搞不懂了。。
趁我现在脑子清新,概念还比较清晰,这里整理一下,也方便将来我又分不清东南西北的时候可以回来复习复习,如果有找不着北的孩子看到这个,希望对你有帮助。
首先我们解释模拟频率,数字频率和采样频率。
假设我们有一小段正弦信号sin(2 π f t),不要考虑任何数字化的东西,这个f就叫做模拟频率。假设f=100,那么我们高中知识知道这个正弦信号的周期是1/f,也就是0.01,也就是说0.01的时间t里面正弦信号变化的一周,我们就知道从0到1里面变化了100次,也就是f次,所以我们才把f叫做频率的。
现在呢,我们要对这个信号进行采样,我们假设每1/1000秒采样一次,也就是说采样周期是0.001秒,那么采样频率就是1000。要信号不失真采样频率是要大于奈奎斯特采样频率的这种事情我就当大家肯定知道好了。
什么是数字频率呢?我们现在假设对这个信号用上面所说的采样频率1000来采样,从t=0采样到t=1,我们知道我们这样就得到了1000个点,假设我们把第一个点称之为a[1],第二个点成为a[2],最后一个点成为a[1000],然后呢,我们知道原本的模拟信号在t=0到1里面变化了100次,由于我们的采样频率大于奈奎斯特频率,所以信号没有失真,所以a[1]到a[1000]里面信号其实也变化了100次,也就是说每10个点相当于原信号变化一个周期。所以说采样之后的数字周期是10,那么数字频率就是1/10。
现在再来看看我们书上强调了很多次的那个公式:
模拟频率=数字频率×采样频率
是不是就吻合了??
上面那个公式换种理解方式就是,原本模拟周期为T模拟的信号经过采样率为f采Hz的采样器后,最后的数字周期就会变成(变大)T数字:
T模拟× f采 = T数字
周期频率换算一下也可以得到上面那个公式。
好了,这三个频率可以理解了,我们就来研究一下DFT(数字傅里叶变换)的频率这个话题。我们知道对于一个信号,如果我们做福利叶变换,那么我们就会得到这个信号的频谱,一个正弦信号,它的福利叶变换,我们知道在它频谱幅度会在正弦信号的频率有一个尖峰(频率的相反数处也有)。但是这个频率是连续的,我们一般做的是数字信号处理,用的是DFT,得到的“傅里叶变换”是一些离散的点,而且点的多少取决于我们想变化的精度。
那么现在假如有一个信号,数字信号,我将其用matlab里的fft做了数字傅里叶变换,得到512个点,那么我怎么知道这些点对应的是什么频率呢??比如说我给的数字信号就一个sin的波形,那么你fft得到的画出来应该是两个尖峰,那么请问这个尖峰对应什么频率呢?
其实回答这个问题不能脱离采样频率这个东西。因为你要相信你给定一个正弦波的数字信号,我肯定可以用随便一个频率的正弦信号,通过精心策划的采样频率得到你那个信号,所以不考虑采样频率的话,是得不到模拟频率这种东西的。
我觉得嘛,理论推导那些东西有兴趣的可以自己去搞一搞,我博客这里写公式不方便,为了更加简洁明了,我直接上代码和图,细细读一下代码,这个问题答案就直接出来了:
Code:
fs = 1000; %采样频率 f = 100; %模拟频率 td = fs/f %数字周期 draw_num = 100;%时域绘画周期数 t = 0:1/fs:draw_num/f; y = sin(2*pi*f*t); X = fftshift(abs(fft(y,1023))); w = (-0.5:1/(length(X)-1):0.5)*fs; size(w) size(X) subplot(211) plot(t,y); subplot(212) plot(w,X)
我们模拟频率f = 100,采样频率fs = 1000,一次构造了一个持续了100个周期的正弦信号y(我们matlab仿真是没有办法做到正弦信号t从-∞到+∞的,向我们这里画出来的其实就是一个理想的正弦信号乘以一个截断的窗函数,反映在频谱上就是窗函数的FFT卷积正弦信号,这样就造成了失真,为了避免这个,我们只能尽可能的使正弦信号持续周期尽可能成)。然后X = fftshift(abs(fft(y,1024)))得到1024个点的数字福利叶变换(fftshift的作用如下:我们知道数字信号的傅里叶变换一定是以2π为周期的,如果直接从fft得到的是0到2π的频谱,移一下之后就可以得到-π到π的了。),频率定标的函数如下面那行代码所示,也就是说我们得到的傅里叶变化的频率实际上是从-fs/2到fs/2的。不管你DFT得到的是多少个点,你的分部一定是在-fs/2到fs/2。再换句话说,你得到的傅里叶变换X,它的数值是X[1],x[2],x[3]…x[1024],那么X[1]对应的频率就是-fs/2,X[1024]对应的就是fs/2。 在换种角度理解这个结果的科学性。我们知道要是信号不失真,采样频率一定要大于信号最高频率的2倍,我们正弦信号的最高频率就是他自己的频率,假设是f,如果我们用正好2f的频率来采样,也就是说fs = 2f,我们再想象一下上面那个图会变成什么样子?两个尖峰就出现在最两边了!!我们看到的只是一个周期,也就是说临界 采样频率下采样的信号,频谱的有效部分正好和另一个周期接壤了,如果我们正好以2f采样,信号正好不可分!!从采样频率也可以理解这个事实了么?? 接下来,再看看数字滤波器,先看下面这两行代码:
fil = dbwavf('db5'); freqz(fil);
功能就是生成一个滤波器的冲击响应,然后看它的频率响应,结果如下:
这是啥子意思捏,一开始学的人都会蒙,就一个0到1或者说0到π的滤波器,它究竟率的是什么频率的呢?其实这个问题,又要回到采样频率和模拟频率那些勾当里面去。。
还是直接上代码和图:
Code:
fs = 300; %采样频率 f = 85; %模拟频率 % td = fs/f; %数字周期 draw_num = 100;%时域绘画周期数 t = 0:1/fs:draw_num/f; y = sin(2*pi*f*t); X = fftshift(abs(fft(y,1023))); w = (-0.5:1/(length(X)-1):0.5)*fs; fil = dbwavf('db5');%随便一个滤波器的冲击响应 yy = conv(fil,y);%滤波 yy = wkeep(yy,length(t));%保持长度 XX = fftshift(abs(fft(yy,1023)));%滤波后频谱 max(XX)/max(X) %查看最高值变化比 应大约等于滤波器频谱该频率下的值 [h wi] = freqz(fil,1,512,fs); w_new = [(-1)*wi(end:-1:2)' wi']'; subplot(321),plot(t,y);title('原信号') subplot(322),plot(w,X);title('原信号频谱') subplot(323),plot(t,yy);title('信号滤波后') subplot(324),plot(w,XX);title('信号滤波后频谱') subplot(325);plot(wi,abs(h));title('滤波器幅度响应') subplot(326);plot(w_new,abs([h(end:-1:2)' h']'));title('滤波器频谱')
代码表示,我们有一个模拟频率f=85的正弦信号,用采样频率fs=300来采样得到一个数字信号,然后DFT算出这个数字信号的频谱,然后我们依旧用上文提到的那个滤波器来研究,那个滤波器和我们的时域信号卷积后,再画出这个滤波后的信号的频谱。那个滤波器的幅度响应如上图左下所示,他的频谱如右下所示,这时我们就看到其实右下这个频谱对应于什么频率是要考虑上采样频率的,我们的采样频率fs=300,所以这个对应的是-150到150。我们的原信号是个单频信号,所以我们要看滤波器与频率的关系很容易。首先我们看左上那个是原正弦信号,滤波后我们看幅度大概只有一半了。我们的频率f是85,我们可以看一下右下角那个图85对应的幅值是多少,可以告诉你,大概就是0.5左右。我们可以再看看正弦信号滤波前后的DFT谱,最高点的幅度也变成了一半左右。
现在终于知道了,我们一个数字信号通过一个滤波器,滤波器响应中1对应的是我们采样频率的一半,0对应的就是0。所以一切的一切其实都是跟采样这种东西有关。单单讨论数字信号不讨论它是怎么采样来的是没有意义的!!
好吧,先到这吧~
【完】
本文内容遵从CC版权协议,转载请注明出自http://www.kylen314.com