MATLAB 2024a 频域分析实战：周期/非周期/系统信号的傅里叶变换与频谱图绘制

频域分析是信号处理领域的核心技能之一，它让我们能够从频率的角度观察和理解信号特性。MATLAB作为工程计算和信号处理的利器，提供了强大的频域分析工具链。本文将带你深入实战，掌握三类典型信号（周期信号、非周期信号和系统响应）的傅里叶分析方法，并完成从理论到代码的完整实现。

1. 周期信号的傅里叶级数展开与频谱分析

周期信号是信号处理中最基础的一类信号，其频域特性可以通过傅里叶级数完美描述。让我们从一个周期为2π的锯齿波信号开始，演示完整的分析流程。

1.1 锯齿波信号的傅里叶系数计算

在MATLAB中，我们可以使用符号计算工具箱来精确求解傅里叶系数。以下是计算过程的关键步骤：

syms t n T; f1 = 4/T * t + 1; % 正半周期表达式 f2 = -4/T * t + 1; % 负半周期表达式 % 计算直流分量a0 a0 = 1/T * (int(f1, -T/2, 0) + int(f2, 0, T/2)); % 计算余弦项系数an ya1 = f1 * cos(n * 2*pi/T * t); ya2 = f2 * cos(n * 2*pi/T * t); an = 2/T * (int(ya1, -T/2, 0) + int(ya2, 0, T/2)); % 计算正弦项系数bn yb1 = f1 * sin(n * 2*pi/T * t); yb2 = f2 * sin(n * 2*pi/T * t); bn = 2/T * (int(yb1, -T/2, 0) + int(yb2, 0, T/2));

注意：对于偶函数，bn理论上应为0，但为了完整性我们仍然计算它。在实际应用中，了解信号的对称性可以显著简化计算。

1.2 频谱图绘制与谐波分析

得到傅里叶系数后，我们可以绘制信号的频谱图，直观展示各频率分量的强度：

n = 1:10; % 谐波次数范围 an1 = subs(an, [T, n], [2*pi, n]); % 代入具体参数值 an2 = [double(a0), double(an1)]; % 包含直流分量 figure; stem([0, n], an2, 'filled'); xlabel('\omega (rad/s)'); ylabel('Amplitude'); title('锯齿波信号的频谱图'); grid on;

通过频谱图可以观察到：

• 锯齿波包含所有整数次谐波
• 谐波幅度随频率增加而衰减
• 偶次谐波分量为零（特定对称性导致）

典型周期信号的频谱特征对比

2. 非周期信号的傅里叶变换实现

非周期信号的频域分析需要用到傅里叶变换。MATLAB提供了fourier()和ifourier()函数进行符号计算，同时也支持数值计算。

2.1 指数衰减信号的频域分析

以单边指数衰减信号f(t)=e^(-3t)u(t)为例：

syms t w; f = exp(-3*t) * heaviside(t); % 定义信号 Fw = fourier(f); % 符号傅里叶变换 % 绘制时域波形和频域特性 figure; subplot(3,1,1); ezplot(f, [0, 2]); title('时域波形'); xlabel('t'); ylabel('f(t)'); % 绘制幅度谱 w_range = -10:0.1:10; Fw_num = subs(Fw, w, w_range); subplot(3,1,2); plot(w_range, abs(Fw_num)); title('幅度谱'); xlabel('\omega'); ylabel('|F(j\omega)|'); % 绘制相位谱 subplot(3,1,3); plot(w_range, angle(Fw_num)); title('相位谱'); xlabel('\omega'); ylabel('\phi(\omega)');

2.2 Sa函数的傅里叶变换对

Sa函数（抽样函数）在信号处理中极为重要，其傅里叶变换是矩形脉冲：

syms t w tao; tao = 3; % 设置参数 Fw = tao * sin(tao*w/2) / (tao*w/2); % 频域表达式 ft = ifourier(Fw, t); % 时域表达式 % 绘制完整的时频分析图 figure; subplot(2,2,1); ezplot(ft, [-5, 5]); title('时域波形'); subplot(2,2,2); ezplot(Fw, [-10, 10]); title('频谱图'); subplot(2,2,3); ezplot(abs(Fw), [-10, 10]); title('幅度谱'); subplot(2,2,4); ezplot(angle(Fw), [-10, 10]); title('相位谱');

常见非周期信号的傅里叶变换对

1. 矩形脉冲 ↔ Sa函数
2. 高斯脉冲 ↔ 高斯函数
3. 单边指数衰减 ↔ 洛伦兹线型
4. 冲激函数 ↔ 常数谱

3. 系统频率响应与失真分析

线性时不变系统(LTI)的频域特性由频率响应函数H(jω)完全描述。我们通过一个具体系统来分析其特性。

3.1 频率响应特性分析

给定系统传递函数H(jω)=(1-jω)/(1+jω)，分析其频率特性：

syms w t; Hw = (1 - 1j*w) / (1 + 1j*w); % 定义系统函数 % 绘制幅频和相频特性 w_range = -5:0.01:5; Hw_num = subs(Hw, w, w_range); figure; subplot(2,1,1); plot(w_range, abs(Hw_num)); title('幅频特性'); xlabel('\omega'); ylabel('|H(j\omega)|'); subplot(2,1,2); plot(w_range, angle(Hw_num)); title('相频特性'); xlabel('\omega'); ylabel('\phi(\omega)');

提示：无失真传输系统需要满足两个条件：1) 幅频响应为常数；2) 相频响应是ω的线性函数。本系统显然不满足第二个条件。

3.2 冲激响应与零状态响应

系统的时域特性可以通过冲激响应来表征：

ht = ifourier(Hw, t); % 计算冲激响应 figure; ezplot(ht, [-2, 5]); title('冲激响应');

当输入为f(t)=e^(-2t)u(t)时，计算系统的零状态响应：

ft = exp(-2*t) * heaviside(t); Fw = fourier(ft); Yw = Fw * Hw; % 频域相乘 yt = ifourier(Yw, t); % 时域响应 figure; subplot(1,2,1); ezplot(ft, [0, 3]); title('输入信号'); subplot(1,2,2); ezplot(yt, [0, 5]); title('系统响应');

4. 综合实战：完整MATLAB脚本实现

将上述分析整合为一个完整的MATLAB脚本，便于实际应用：

%% 周期信号分析 - 锯齿波 T = 2*pi; syms t n; f1 = 4/T*t + 1; f2 = -4/T*t + 1; % 计算傅里叶系数 a0 = 1/T*(int(f1,t,-T/2,0)+int(f2,t,0,T/2)); an = 2/T*(int(f1*cos(n*2*pi/T*t),t,-T/2,0)+int(f2*cos(n*2*pi/T*t),t,0,T/2)); bn = 2/T*(int(f1*sin(n*2*pi/T*t),t,-T/2,0)+int(f2*sin(n*2*pi/T*t),t,0,T/2)); % 绘制频谱 n_harmonics = 1:10; an_vals = double(subs(an, n, n_harmonics)); figure; stem([0, n_harmonics], [double(a0), an_vals], 'filled'); title('锯齿波频谱'); xlabel('谐波次数'); ylabel('幅度'); %% 非周期信号分析 - 指数衰减信号 f = exp(-3*t)*heaviside(t); Fw = fourier(f); w = -10:0.1:10; Fw_num = subs(Fw, w, w); figure; subplot(2,1,1); plot(w, abs(Fw_num)); title('幅度谱'); subplot(2,1,2); plot(w, angle(Fw_num)); title('相位谱'); %% 系统分析 Hw = (1-1j*w)/(1+1j*w); ht = ifourier(Hw, t); % 计算特定输入的响应 ft = exp(-2*t)*heaviside(t); yt = ifourier(fourier(ft)*Hw, t);

工程应用中的实用技巧

1. 对于数值计算，优先使用fft而非符号计算，效率更高
2. 分析长信号时，采用分段加窗的方法减少频谱泄漏
3. 系统辨识中，可以通过频响函数估计系统参数
4. 使用fftshift函数将零频分量移到频谱中心
5. 对于实时处理，考虑使用Goertzel算法计算特定频率分量