用MATLAB的FFT函数对已测时域数据进行傅立叶变换

用MATLAB的FFT函数对已测时域数据进行傅立叶变换

一、原理

FFT：快速傅立叶变换
1.相关代码用法

Y = fft(X) 
%计算X的离散傅立叶变换，实质是返回一个向量或多个向量的傅立叶变换Y = fft(X,n) 
%n指定Y的长度，如果未指定，默认X的长度与Y的长度相等Y = fft(X,n,dim) 
%如果X是矩阵，则fft(X,n,2)返回每行的n点傅里叶变换

2.示例代码的理解

比较时域和频域中的余弦波。
指定信号的参数，采样频率为 1kHz，信号持续时间为 1 秒。

Fs = 1000;                % Sampling frequency
T = 1/Fs;                 % Sampling period
L = 1000;                 % Length of signal
t = (0:L-1)*T;            % Time vector

注释：L代表采样个数，t是一个离散向量，代表每隔采样时间T进行一次采样，总共采L次样，所以总的时间区间为(0:L-1)*T。

创建一个矩阵，其中每一行代表一个频率经过缩放的余弦波。结果 X 为 3×1000 矩阵。第一行的波频为 50，第二行的波频为 150，第三行的波频为 300。

x1 = cos(2*pi*50*t);          % First row wave
x2 = cos(2*pi*150*t);         % Second row wave
x3 = cos(2*pi*300*t);         % Third row wave
X = [x1; x2; x3];

x1，x2，x3是由离散序列t产生的离散向量，X是三个行向量组成的矩阵，其列的长度取决于定义的L的长度。

在单个图窗中按顺序绘制 X 的每行的前 100 个项，并比较其频率。

for i = 1:3subplot(3,1,i) %把画图空间分成三行一列，在第i行进行画图plot(t(1:100),X(i,1:100))title(['Row ',num2str(i),' in the Time Domain'])
end

出于算法性能的考虑，fft 允许用尾随零填充输入。在这种情况下，用零填充 X 的每一行，以使每行的长度为比当前长度大的下一个最小的 2 的次幂值。使用 nextpow2 函数定义新长度。

n = 2^nextpow2(L);
% 举个例子，如果x等于100，则y=7，因为2的7次方==128，而128是所有大于100的，二的整数次幂数字中最小的一个dim = 2;
%指定dim参数沿X的行（即对每个信号）使用fft。dim=2时沿行返回每个离散序列的傅立叶变换Y = fft(X,n,dim);P2 = abs(Y/L); 
% 计算双侧频谱 P2
P1 = P2(:,1:n/2+1);	
% 将P2的前半段信号赋给P1，P1即是我们关心的部分
P1(:,2:end-1) = 2*P1(:,2:end-1);
%计算每个信号的双侧频谱和单侧频谱。for i=1:3subplot(3,1,i)plot(0:(Fs/n):(Fs/2-Fs/n),P1(i,1:n/2))title(['Row ',num2str(i),' in the Frequency Domain'])
end
%在频域内，为单个图窗中的每一行绘制单侧幅值频谱。

3.理解学习
参考资料：Matlab中fft函数的使用与原理

（1）Fn = （n-1）* Fs /N

Fn是第n点所表示的真实频率值。当然，n只取前一半的点就足够了。这样，可以达到的频率分辨率即为Fs/N。

（2）幅值的处理

作FFT分析时，幅值大小与输入点数有关，要得到真实的幅值大小，需要将变换后的结果除以N。且由于零频在双边谱中本没有被一分为二，所以对于零频外的点还有乘以2，得到的才是真实的频率幅值。

（3）基频

​若分析数据时长为T，则分析结果的基频就是f0=1/T，分析结果的频率序列为[0:N-1]*f0

二、运用

从实验室里的NI里导出了一个时域波形的数据，是对正弦信号进行采样保持后从采样保持器里输出的信号，实验目的是对时域的数据进行fft处理的到频谱图，然后与NI频谱分析仪里面的频谱图进行对比，观察是否一致。

待处理数据是时间区间为[24.890153,24.910153]，采样信号频率为1000Hz，被采样信号频率为100Hz的正弦信号。

总的时间区间的长度为 t = 24.91 - 24.89 = 0.02 s
采样频率为1000Hz，采样周期T = 0.001 s
故在这个时间长度里总共能进行L = 0.02/0.001 + 1 = 21次采样
但是总共导出了2500个数据，从这2500个数据中等间隔抽出21个数据，作为傅立叶变换的离散序列。2500/20 = 125。（如果t = 0.2s要考虑：根据抽样保持的特殊性，可以选择每隔12.5*k，然后在excel中截去小数向上取整数 （=CEILING(A1,1））来选取每一个离散点，即=CEILING(A1,1），然后总的操作方法就是在excel表格里面输入：=OFFSET(A$1,CEILING(ROW(A1)*12.5,1),)，然后再让时间区间里面的每个时间点减去起始记录时间时间24.91，得到新的时间区间为[0,0.2] ）

得到第13列为时间离散序列t，14列为对应的采样值。

对应的代码如下

clc;clear;close all;
load('EP01')x1 = EP01(1:21,13);
y1 = EP01(1:21,14);  %需要做傅立叶变换的一组数据
N = 21;   %信号的长度，能不能等价于待处理的数据有多少个
Fs = 1000;yk=fft(y1);      % 傅立叶变换
P2 = abs(yk/N); 
% 计算双侧频谱 P2
P1 = P2(1:N/2+1);	
% 将P2的前半段信号赋给P1，P1即是我们关心的部分
P1(2:end-1) = 2*P1(2:end-1);
%计算每个信号的双侧频谱和单侧频谱。
plot(0:round(Fs/N):round(Fs/2-Fs/N),P1(1:N/2))title('FFT变换')
legend('matlab进行输出信号的fft变换')
xlabel('频率/Hz')
ylabel('归一化幅度')

得到的频谱图如下

如果再经过一个低通滤波器，应该会有更完美的频谱特性。

接下来对采样保持输出器的信号做频谱分析。经过采样保持器后的波形如下：

同理：
NI的采样频率为Ts = 125000
N = 2500
代入代码即可

clc;clear;close all;
load('EP01')
x = EP01(1:400,11);
y = EP01(1:400,12);
% y = EP01(1:2500,9);x1 = EP01(1:2500,10);
y1 = EP01(1:2500,9);  %需要做傅立叶变换的一组数据
N = 2500;   %信号的长度，能不能等价于待处理的数据有多少个
Fs = 125000;subplot(211);   %分为两行一列在第一个位置画图
plot(x,y)
% xlim([24.89,24.91])
grid on 
legend('采样保持器的输出信号的频域波形')
xlabel('频率/Hz')
ylabel('dB')yk=fft(y1);      % 傅立叶变换
P2 = abs(yk/N); 
% 计算双侧频谱 P2
P1 = P2(1:N/2+1);	
% 将P2的前半段信号赋给P1，P1即是我们关心的部分
P1(2:end-1) = 2*P1(2:end-1);
%计算每个信号的双侧频谱和单侧频谱。
subplot(212);plot(0:round(Fs/N):round(Fs/2-Fs/N),abs(P1(1:N/2)))
% plot(x1,abs(yk));
xlim([0;40000])
title('FFT变换')
legend('matlab进行输出信号的fft变换')
xlabel('频率/Hz')
ylabel('归一化幅度')

画出来的频谱图和实际导出来的不太一样，找到原因了后续改正。