FFT

来源:互联网 发布:上瘾网络剧北京发布会 编辑:程序博客网 时间:2024/06/16 07:11

快速傅立叶变换(FFT)

函数FFT和IFFT用于快速傅立叶变换和逆变换。下面介绍这些函数。

函数的一种调用格式为         y=fft(x)

其中,x是序列,y是序列的FFT,x可以为一向量或矩阵,若x为一向量,y是x的FFT。且和x相同长度。若x为一矩阵,则y是对矩阵的每一列向量进行FFT。

如果x长度是2的幂次方,函数fft执行高速基-2FFT算法;否则fft执行一种混合基的离散傅立叶变换算法,计算速度较慢。

函数FFT的另一种调用格式为          y=fft(x,N)

式中,x,y意义同前,N为正整数。

函数执行N点的FFT。若x为向量且长度小于N,则函数将x补零至长度N。若向量x的长度大于N,则函数截短x使之长度为N。若x 为矩阵,按相同方法对x进行处理。

经函数fft求得的序列y一般是复序列,通常要求其幅值和相位。MATLAB提供求复数的幅值和相位函数:abs,angle,这些函数一般和FFT同时使用。

函数abs(x)用于计算复向量x的幅值,函数angle(x)用于计算复向量的相角,介于 和 之间,以弧度表示。

函数unwrap(p)用于展开弧度相位角p ,当相位角绝对变化超过 时,函数把它扩展至 。

实验1.

[plain] view plain copy
  1. % x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t)。采样频率fs=100Hz,分别绘制N=128、1024点幅频图  
  2. clf;  
  3. fs=100;N=128;   %采样频率和数据点数  
  4. n=0:N-1;t=n/fs;   %时间序列  
  5. x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t); %信号  
  6. y=fft(x,N);    %对信号进行快速Fourier变换  
  7. mag=abs(y);     %求得Fourier变换后的振幅  
  8. f=n*fs/N;    %频率序列  
  9. subplot(321),plot(x);  axis([0 N-1 min(x) max(x) ])   
  10. xlabel('时间/t');ylabel('振幅');title('信号');grid on;  
  11. subplot(322),plot(y);   
  12. xlabel('频率/Hz');ylabel('振幅');title('做FFT后');grid on;  
  13. subplot(323),plot(f,mag);   %绘出随频率变化的振幅  
  14. xlabel('频率/Hz');ylabel('振幅');title('N=128');grid on;  
  15. subplot(324),plot(f(1:N/2),mag(1:N/2)); %绘出Nyquist频率之前随频率变化的振幅  
  16. xlabel('频率/Hz');ylabel('振幅');title('N=128');grid on;  
  17. %对信号采样数据为1024点的处理  
  18. fs=100;N=1024;n=0:N-1;t=n/fs;  
  19. x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t); %信号  
  20. y=fft(x,N);   %对信号进行快速Fourier变换  
  21. mag=abs(y);   %求取Fourier变换的振幅  
  22. f=n*fs/N;  
  23. subplot(325),plot(f,mag); %绘出随频率变化的振幅  
  24. xlabel('频率/Hz');ylabel('振幅');title('N=1024');grid on;  
  25. subplot(326)  
  26. plot(f(1:N/2),mag(1:N/2)); %绘出Nyquist频率之前随频率变化的振幅  
  27. xlabel('频率/Hz');ylabel('振幅');title('N=1024');grid on;  

所得结果为:

从图中可以看出:x信号共有15HZ和40HZ两个,在做FFT后都表现了出来,并且它们是以f=50HZ对称的。

振幅的大小与所用采样点数有关,采用128点和1024点的相同频率的振幅是有不同的表现值,但在同一幅图中,40Hz与15Hz振动幅值之比均为4:1,与真实振幅0.5:2是一致的。为了与真实振幅对应,需要将变换后结果乘以2除以N。

实验2.

[plain] view plain copy
  1. Fs = 1000;                    % 采样频率  
  2. T = 1/Fs;                     % 采样时间  
  3. L = 1000;                     % 总的采样点数  
  4. t = (0:L-1)*T;                % 时间序列(时间轴)  
  5. %产生一个幅值为0.7频率为50HZ正弦+另外一个信号的幅值为1频率为120Hz的正弦信号  
  6. x = 0.7*sin(2*pi*50*t) + sin(2*pi*120*t);   
  7. y = x + 2*randn(size(t));     % 混入噪声信号  
  8. subplot(211)  
  9. plot(Fs*t(1:50),y(1:50))      %画出前50个点  
  10. title('Signal Corrupted with Zero-Mean Random Noise')  
  11. xlabel('time (milliseconds)')  
  12. NFFT = 2^nextpow2(L); % 求得最接近总采样点的2^n,这里应该是2^10=1024  
  13. Y = fft(y,NFFT)/L;    %进行fft变换(除以总采样点数,是为了后面精确看出原始信号幅值)  
  14. f = Fs/2*linspace(0,1,NFFT/2+1);%频率轴(只画到Fs/2即可,由于y为实数,后面一半是对称的)  
  15. % 画出频率幅度图形,可以看出50Hz幅值大概0.7,120Hz幅值大概为1.  
  16. subplot(212)  
  17. plot(f,2*abs(Y(1:NFFT/2+1)))   
  18. title('Single-Sided Amplitude Spectrum of y(t)')  
  19. xlabel('Frequency (Hz)')  
  20. ylabel('|Y(f)|')  


实验3.

对语音信号做FFT,

[plain] view plain copy
  1. [x,fs,bits]=wavread('test1.wav');%读入语音信号  
  2. sound(x,fs,bits); %回放语音  
  3. X=fft(x,1024);  
  4. magX=abs(X);  
  5. angX=angle(X);  
  6. subplot(221);plot(x);title('原始信号波形');  
  7. subplot(222);plot(X); title('原始信号频谱');  
  8. subplot(223);plot(magX);title('原始信号幅值');  
  9. subplot(224);plot(angX);title('原始信号相位');  

结果: