FFT
来源：互联网 发布：上瘾网络剧北京发布会 编辑：程序博客网 时间：2024/06/16 01:37
快速傅立叶变换（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.
- % x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t)。采样频率fs=100Hz，分别绘制N=128、1024点幅频图
- clf;
- fs=100;N=128; %采样频率和数据点数
- n=0:N-1;t=n/fs; %时间序列
- x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t); %信号
- y=fft(x,N); %对信号进行快速Fourier变换
- mag=abs(y); %求得Fourier变换后的振幅
- f=n*fs/N; %频率序列
- subplot(321),plot(x); axis([0 N-1 min(x) max(x) ])
- xlabel('时间/t');ylabel('振幅');title('信号');grid on;
- subplot(322),plot(y);
- xlabel('频率/Hz');ylabel('振幅');title('做FFT后');grid on;
- subplot(323),plot(f,mag); %绘出随频率变化的振幅
- xlabel('频率/Hz');ylabel('振幅');title('N=128');grid on;
- subplot(324),plot(f(1:N/2),mag(1:N/2)); %绘出Nyquist频率之前随频率变化的振幅
- xlabel('频率/Hz');ylabel('振幅');title('N=128');grid on;
- %对信号采样数据为1024点的处理
- fs=100;N=1024;n=0:N-1;t=n/fs;
- x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t); %信号
- y=fft(x,N); %对信号进行快速Fourier变换
- mag=abs(y); %求取Fourier变换的振幅
- f=n*fs/N;
- subplot(325),plot(f,mag); %绘出随频率变化的振幅
- xlabel('频率/Hz');ylabel('振幅');title('N=1024');grid on;
- subplot(326)
- plot(f(1:N/2),mag(1:N/2)); %绘出Nyquist频率之前随频率变化的振幅
- 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.
- Fs = 1000; % 采样频率
- T = 1/Fs; % 采样时间
- L = 1000; % 总的采样点数
- t = (0:L-1)*T; % 时间序列（时间轴）
- %产生一个幅值为0.7频率为50HZ正弦+另外一个信号的幅值为1频率为120Hz的正弦信号
- x = 0.7*sin(2*pi*50*t) + sin(2*pi*120*t);
- y = x + 2*randn(size(t)); % 混入噪声信号
- subplot(211)
- plot(Fs*t(1:50),y(1:50)) %画出前50个点
- title('Signal Corrupted with Zero-Mean Random Noise')
- xlabel('time (milliseconds)')
- NFFT = 2^nextpow2(L); % 求得最接近总采样点的2^n,这里应该是2^10=1024
- Y = fft(y,NFFT)/L; %进行fft变换（除以总采样点数，是为了后面精确看出原始信号幅值）
- f = Fs/2*linspace(0,1,NFFT/2+1);%频率轴(只画到Fs/2即可，由于y为实数，后面一半是对称的)
- % 画出频率幅度图形，可以看出50Hz幅值大概0.7,120Hz幅值大概为1.
- subplot(212)
- plot(f,2*abs(Y(1:NFFT/2+1)))
- title('Single-Sided Amplitude Spectrum of y(t)')
- xlabel('Frequency (Hz)')
- ylabel('|Y(f)|')
实验3.
对语音信号做FFT，
- [x,fs,bits]=wavread('test1.wav');%读入语音信号
- sound(x,fs,bits); %回放语音
- X=fft(x,1024);
- magX=abs(X);
- angX=angle(X);
- subplot(221);plot(x);title('原始信号波形');
- subplot(222);plot(X); title('原始信号频谱');
- subplot(223);plot(magX);title('原始信号幅值');
- subplot(224);plot(angX);title('原始信号相位');
