一文搞懂Matlab时频傅里叶转换——修正后重发

昨天发表的《一文搞懂Matlab时频傅里叶转换》一文中存在严重的理论错误，感谢CAE中学生的热心群友提醒，错误修正的地方将在下文详细说明。

众所周知，振动中会用到响应谱，响应谱分析的主要用途是代替瞬态动力学的时程分析，如下图，将一个复杂的时域信号使用拟合的方法，分离成众多的谐波信号，即正弦或余弦信号，由于频谱图省略了各谐波的相位信息，只记录它们的频率与幅值，所以是近似且快速的计算手段。这个内容在教主的动力学课程中说得很详细了。

严格意义上，响应谱分析中的边界条件是一种响应，因为输出的反应谱是位移（速度、加速度或力）与频率之间的关系，这种关系称为频谱。

但是大多数时候，我们要计算的频谱并不是标准手册上有的，而是根据我们自身产品在冲击或振动试验中测得的。一般使用传感器测试产品某处的时间历程响应，即时间与位移（或速度、加速度）的时域信号。如果是在固定约束处，也可直接使用振动台的时域信号，而无需传感器。将这个时域信号转换为频域信号，导入有限元软件中对产品进行仿真计算。如果振动台控制器自带频谱输出，那么都不需要时频转换了。

工作中，甲方也不会给你配一个人专门进行信号处理，往往甲方自己也不清楚边界或加载情况，更有甚者试验都让仿真者参与设计，所以不要想着有人给你个频谱，你拿进ansys/abaqus中咔咔一顿操作交差完事，这么轻松为啥甲方自己不搞呢。

在有限元教材中不会讲Matlab的时频转换，在Matlab的教材中又不会讲振动，这个最重要的环节似乎被各位大佬故意遗忘了，我们本文中要探讨的正是这个中间被忽略的部分。

上面划线这部分存在理论错误。先看周老师的说明：

从严格意义上来说，响应谱分析中的边界条件不是传统意义上的载荷，而是一种响应，因为输入的反应谱是位移、速度、加速度或力与频率之间的关系，这种关系被称为频谱。例如:在振动试验台上安装4个单自由度的弹簧质量系统,频率分别为f1、f2、f3、f4，且f1＜f2＜f3＜f4。当振动台以频率f1激振时,可得4个系统的位移响应图，再叠加f3频率激振效果，可得4个系统的位移响应图，最后叠加所有频率的综合激振效果，得到的曲线为频谱，如图 4-1-2 所示。  

个人理解，振动台只是起到扫频作用，比如从0~1000HZ扫频，在某些频率下，产品会发生共振而导致振幅（位移、速度或加速度）变大，记录下0~1000HZ下产品所有的振幅，形成响应谱图。  

所以可见，响应谱是产品在连续频率的输入载荷下表现出的自身响应特性。同样0~1000HZ扫频，产品A与产品B的响应谱可能完全不同。

本文使用的是16版本的Matlab，按照Matlab桀骜不驯的性格，每年升级都会悄悄咪 咪地改动几个命令，所以不能保证适用于其他版本的Matlab。

Step1 数据分析    

如图，打开ZXS.csv文件，可以看到这个信号共1s，每个采样数据间隔0.01s，即采样率100。这里一共有100个数据，注意数据的数量必须是偶数个，数据最上方不能有表头，否则需要修改后文的代码。

Step2 Matlab转换。

首先将工作目录设置到ZXS.csv所在的文件夹。

输入如下代码，%之后的文字为注释。

需要注意的是，需要核对输入文件名称ZXS.csv是否与我们保存的数据文件名相同，核对采样率是否正确。

clc;clear all;close all %清理屏幕与缓存

data = xlsread('ZXS.csv'); %19版本读取命令为data = readmatrix('ZXS.csv')    

t = data(:,1);      %对应第一列

y1 = data(:,2);     %对应第二列

figure

y=y1;         % 读取时域数据

Fs=100;     % 采集频率，1秒采集了多少个点（需要自行根据数据更改）

T=1/Fs;       % 采集时间间隔

N=length(y);  % 采集信号的长度

t=(0:1:N-1)*T;   % 定义整个采集时间点

t=t';            % 转置成列向量

% fft变换

Y=fft(y);          % Y为fft变换结果，复数向量

Y=Y(1:N/2+1);      % 只看变换结果的一半即可

A=abs(Y);          % 复数的幅值（模）

f=(0:1:N/2)*Fs/N;  % 生成频率范围

f=f';              % 转置成列向量

% 幅值修正

A_adj=zeros(N/2+1,1);

A_adj(1)=A(1)/N;      % 频率为0的位置

A_adj(end)=A(end)/N;   % 频率为Fs/2的位置

A_adj(2:end-1)=2*A(2:end-1)/N;

%写出数据    

dataa(:,1)=f;               %频率

dataa(:,2)=A_adj;           %幅值

str = ['DXS.csv'];     %写出文件名称

xlswrite(str,dataa);    %19版本为  writematrix(str,dataa)

% 绘制频率幅值图和频谱图

subplot(2,1,1)

plot(t,y,'LineWidth',1)

xlabel('时间')

ylabel('信号值')

title('时域信号')

subplot(2,1,2)

plot(f,A_adj,'LineWidth',3)

xlabel('频率(Hz)')

ylabel('幅值(修正后)')

title('FFT变换幅值图')

输入后，点击上方运行工具，将自动导出转换的频域信号DXS.csv，并生成图像，上面显示的是输入的时域图，下面显示的是输出的频域图。可以看到，这个型号是由三个信号组成，分别是频率=1幅值=5的信号1，频率=8幅值=1的信号2，频率=10幅值=0.5的信号3。   

到输出的DXS.csv文件中，插入散点图能更清晰看到组成。

可以看到，频域信号底部前后各扩展了1hz，形成了三角形，这是因为频谱分辨率导致，分辨率=1/总时间，我们的信号总时间是1s，所以分辨率是1Hz。

现在，你可以将自己的时频信号导入测试了。

动力学的仿真实操其实很简单，与静力学无异，难就难在理论的解读，我也感觉到自身理论的短板，也请大家多大指教批评。

本文特别感谢卧虎藏龙的CAE中学生的初一（2）班的jiff、黄荣等大佬的热心解答。

参考书籍：《ANSYS  Workbench  有限元分析实例详解  动力学》—周炬、苏金英

参考视频：《ANSYS高级视频教程》——张晔