MATLAB快速傅里叶变换(fft)函数详解
MATLAB快速傅里叶变换(fft)函数详解
调用:
1. Y = fft(y);
- Y = fft(y,N);
式中,y是序列,Y是序列的快速傅里叶变换。y可以是一向量或矩阵,若y为向量,则Y是y的FFT,并且与y具有相同的长度。若y为一矩阵,则Y是对矩阵的每一列向量进行FFT。
说明:
函数fft返回值的数据结构具有对称性
根据采样定理,fft能分辨的最高频率为采样频率的一半(即Nyquist频率),函数fft返回值是以Nyqusit频率为轴对称的,Y的前一半与后一半是复数共轭关系。
幅值
作FFT分析时,幅值大小与输入点数有关,要得到真实的幅值大小,只要将变换后的结果乘以2除以N即可(但此时零频—直流分量—的幅值为实际值的2倍)。对此的解释是:Y除以N得到双边谱,再乘以2得到单边谱(零频在双边谱中本没有被一分为二,而转化为单边谱过程中所有幅值均乘以2,所以零频被放大了)。
基频
若分析数据时长为T,则分析结果的基频就是f0=1/T,分析结果的频率序列为[0:N-1]*f0
执行N点FFT
在调用格式2中,函数执行N点FFT。若y为向量且长度小于N,则函数将y补零至长度N,若向量y的长度大于N,则函数截断y使之长度为N。
注意:
使用N点FFT时,若N大于向量y的长度,将给频谱分析结果带来变化,应该特别注意。
结论:
使用N点FFT时,不应使N大于y向量的长度,否则将导致频谱失真。
例子程序:
clear all %清除内存所有变量
close all %关闭所有打开的图形窗口
%% 执行FFT点数与原信号长度相等(100点)
% 构建原信号
N=100; % 信号长度(变量@@@@@@@)
Fs=1; % 采样频率
dt=1/Fs; % 采样间隔
t=[0:N-1]*dt; % 时间序列
xn=cos(2*pi*0.24*[0:99])+cos(2*pi*0.26*[0:99]);
xn=[xn,zeros(1,N-100)]; % 原始信号的值序列
subplot(3,2,1) % 变量@@@@@@@
plot(t,xn) % 绘出原始信号
xlabel('时间/s'),title('原始信号(向量长度为100)') % 变量@@@@@@@
% FFT分析
NN=N; % 执行100点FFT
XN=fft(xn,NN)/NN; % 共轭复数,具有对称性
f0=1/(dt*NN); % 基频
f=[0:ceil((NN-1)/2)]*f0; % 频率序列
A=abs(XN); % 幅值序列
subplot(3,2,2),stem(f,2*A(1:ceil((NN-1)/2)+1)),xlabel('频率/Hz') % 绘制频谱(变量@@@@@@@)
axis([0 0.5 0 1.2]) % 调整坐标范围
title('执行点数等于信号长度(单边谱100执行点)'); % 变量@@@@@@@
%% 执行FFT点数大于原信号长度
% 构建原信号
N=100; % 信号长度(变量@@@@@@@)
Fs=1; % 采样频率
dt=1/Fs; % 采样间隔
t=[0:N-1]*dt; % 时间序列
xn=cos(2*pi*0.24*[0:99])+cos(2*pi*0.26*[0:99]);
xn=[xn,zeros(1,N-100)]; % 原始信号的值序列
subplot(3,2,3) % 变量@@@@@@@
plot(t,xn) % 绘出原始信号
xlabel('时间/s'),title('原始信号(向量长度为100)') % 变量@@@@@@@
% FFT分析
NN=120; % 执行120点FFT(变量@@@@@@@)
XN=fft(xn,NN)/NN; % 共轭复数,具有对称性
f0=1/(dt*NN); % 基频
f=[0:ceil((NN-1)/2)]*f0; % 频率序列
A=abs(XN); % 幅值序列
subplot(3,2,4),stem(f,2*A(1:ceil((NN-1)/2)+1)),xlabel('频率/Hz') % 绘制频谱(变量@@@@@@@)
axis([0 0.5 0 1.2]) % 调整坐标范围
title('执行点数大于信号长度(单边谱120执行点)'); % 变量@@@@@@@
%% 执行FFT点数与原信号长度相等(120点)
% 构建原信号
N=120; % 信号长度(变量@@@@@@@)
Fs=1; % 采样频率
dt=1/Fs; % 采样间隔
t=[0:N-1]*dt; % 时间序列
xn=cos(2*pi*0.24*[0:99])+cos(2*pi*0.26*[0:99]);
xn=[xn,zeros(1,N-100)]; % 原始信号的值序列
subplot(3,2,5) % 变量@@@@@@@
plot(t,xn) % 绘出原始信号
xlabel('时间/s'),title('原始信号(向量长度为120)') % 变量@@@@@@@
% FFT分析
NN=120; % 执行120点FFT(变量@@@@@@@)
XN=fft(xn,NN)/NN; % 共轭复数,具有对称性
f0=1/(dt*NN); % 基频
f=[0:ceil((NN-1)/2)]*f0; % 频率序列
A=abs(XN); % 幅值序列
subplot(3,2,6),stem(f,2*A(1:ceil((NN-1)/2)+1)),xlabel('频率/Hz') % 绘制频谱(变量@@@@@@@)
axis([0 0.5 0 1.2]) % 调整坐标范围
title('执行点数等于信号长度(单边谱120执行点)'); % 变量@@@@@@@
自相关分析数据处理代码
正弦信号仿真
load('matlab4.mat')
N=1000; fs=100; %数据长度和采样频率
n=0:N-1;t=n/fs; %时间序列
maxlags=500; %延迟样点数
x=2*sin(2*pi*2*t); %原始正弦信号
y=awgn(x,10*log10(0.7));
[c1,lags]=xcorr(x,maxlags,'unbiased');
[c2,lags]=xcorr(y,maxlags,'unbiased'); figure(1)
subplot(2,2,1),plot(t,x); %绘原始信号x
xlabel('时间/s');ylabel('x(t)');title('正弦周期信号');grid on;
subplot(2,2,2),plot(lags/fs,c1);
xlabel('时间/s');ylabel('Rx(t)');title('正弦周期信号自相关');grid on;
subplot(2,2,3),plot(t,y);
xlabel('时间/s');ylabel('y(t)');title('带噪声周期信号');grid on;
subplot(2,2,4),plot(lags/fs,c2);
xlabel('时间/s');ylabel('Ry(t)');title('带噪声周期信号的自相关');grid on
y1=fft(x,N); %对信号进行傅里叶变换
yy1=abs(y1); %求得傅里叶变换后的振幅
yy1=yy1*2/N; %幅值处理
f=n*fs/N; figure(2)
subplot(2,2,1),plot(f(1:N/2),yy1(1:N/2));
xlabel('频率/HZ');ylabel('振幅');title('正弦周期信号的傅里叶变换');grid on;
y2=fft(y,N); %对正弦周期信号进行傅里叶变换
yy2=abs(y2); %求得傅里叶变换后的振幅
yy2=yy2*2/N; %幅值处理
f=n*fs/N; %频率序列
subplot(2,2,2),plot(f(1:N/2),yy2(1:N/2))
xlabel('频率/HZ');ylabel('振幅');title('加噪正弦周期信号傅里叶变换');grid on;
y4=fft(c1,N); %对正弦周期信号进行傅里叶变换
yy4=abs(y4); %求得傅里叶变换后的振幅
yy4=yy4*2/N; %幅值处理
f=n*fs/N; %频率序列
subplot(2,2,3),plot(f(1:N/2),yy4(1:N/2));
xlabel('频率/HZ');ylabel('振幅');title('正弦周期信号自相关傅里叶变换');grid on;
y3=fft(c2,N); %对加噪信号自相关进行傅里叶变换
yy3=abs(y3); %求得傅里叶变换后的振幅
yy3=yy3*2/N; %幅值处理
f=n*fs/N; %频率序列
subplot(2,2,4),plot(f(1:N/2),yy3(1:N/2));
xlabel('频率/HZ');ylabel('振幅');title('加噪正弦周期信号自相关傅里叶变换');grid on;
实测数据的处理
x=torque4;%用M代替扭矩
N=1:1:100001;fs=1000; figure(1)
subplot(2,2,1);plot(N,x);%绘制离散散点图
grid on
xlabel('序列号N');ylabel('序列值M');title('序列值M-序列号N图像');
t=N/fs;%时间序列
subplot(2,2,2);plot(t,x);%绘制随时间的变化图
grid on
xlabel('时间t/s');ylabel('扭矩M/N*m');title('扭矩M-时间t图像');
torque(26690:1:27770)=torque4(26690:1:27770)-mean(torque4(26690:1:27770));
M=torque(26690:1:27770); %截取平稳运行段的扭矩数据并去均值
t1=(26690:1:27770)/fs; %截取后的时间序列
subplot(2,2,3),plot(t1,M);grid on;
xlabel('时间t/s');ylabel('扭矩M/N*m');title('截取平稳扭矩M-时间t1图像');
N1=1081;n=1:N1;Lag=500;
[c,lags]=xcorr(M,Lag,'unbiased');
subplot(2,2,4),plot(lags/fs,c);grid on;
xlabel('时间t/s');ylabel('Rm(t)');title('截取平稳扭矩自相关图像');
y1=fft(M,N1); %对自相关进行傅里叶变换
yy1=abs(y1); %求得傅里叶变换后的振幅
yy1=yy1*2/N1; %幅值处理
f=n*fs/N1; figure(2)
subplot(2,1,1),plot(f(1:N1/2),yy1(1:N1/2));grid on;
xlabel('频率/HZ');ylabel('振幅');title('截取信号的傅里叶变换');
y2=fft(c,N1); %对加噪信号进行傅里叶变换
yy2=abs(y2); %求得傅里叶变换后的振幅
yy2=yy2*2/N1; %幅值处理
f=n*fs/N1; %频率序列
subplot(2,1,2),plot(f(1:N1/2),yy2(1:N1/2))
xlabel('频率/HZ');ylabel('振幅');title('截取平稳运行段的自相关傅里叶变换');grid on;
自振频率的测量采用自相关的方法:
y1=fft(c,N1); %对自相关进行傅里叶变换
yy1=abs(y1); %求得傅里叶变换后的振幅
yy1=yy1*2/N1; %幅值处理
f=n*fs/N1; %频率序列
figure(2)
subplot(2,2,1),plot(f(1:N1/2),yy1(1:N1/2));grid on;
xlabel('频率/HZ');ylabel('振幅');title('截取平稳运行段的自相关傅里叶变换');
y2=fft(M,N1); %对加噪信号进行傅里叶变换
yy2=abs(y2); %求得傅里叶变换后的振幅
yy2=yy2*2/N1; %幅值处理
f=n*fs/N1; %频率序列
subplot(2,2,2),plot(f(1:N1/2),yy2(1:N1/2))
xlabel('频率/HZ');ylabel('振幅');title('截取信号的傅里叶变换');grid on