从MATLAB代码实战看FS、FT、DFS、DTFS、DTFT的区别与应用

从MATLAB代码实战看FS、FT、DFS、DTFS、DTFT的区别与应用在信号处理领域傅里叶分析是一把打开频域大门的金钥匙。但对于许多工程师和学生来说各种傅里叶变换的变体——FS傅里叶级数、FT傅里叶变换、DFS离散傅里叶级数、DTFS离散时间傅里叶级数和DTFT离散时间傅里叶变换——常常让人感到困惑。本文将通过MATLAB代码实例直观展示这些变换在实际应用中的差异和适用场景。1. 理论基础与核心概念傅里叶分析的核心思想是将信号分解为不同频率的正弦波组合。不同类型的信号连续/离散、周期/非周期对应不同的傅里叶分析方法信号类型与变换对应关系信号类型连续信号离散信号周期信号FS (傅里叶级数)DFS/DTFS (离散傅里叶级数)非周期信号FT (傅里叶变换)DTFT (离散时间傅里叶变换)表1不同类型信号对应的傅里叶分析方法1.1 连续信号分析对于连续时间信号我们有两种主要分析方法傅里叶级数(FS)适用于周期连续信号% 生成周期为T的方波信号 T 2*pi; % 周期 t linspace(0, 3*T, 1000); % 时间向量 f square(t); % 方波信号傅里叶变换(FT)适用于非周期连续信号% 生成非周期信号高斯脉冲 t linspace(-5, 5, 1000); f exp(-t.^2); % 高斯脉冲提示MATLAB中的fft函数实际上计算的是DFT离散傅里叶变换要分析连续信号需要适当选择采样率。1.2 离散信号分析对于离散时间信号我们也有对应的分析方法离散傅里叶级数(DFS/DTFS)适用于周期离散信号% 生成周期离散信号 N 20; % 周期 n 0:3*N-1; % 离散时间索引 x mod(n, N) N/2; % 周期方波离散时间傅里叶变换(DTFT)适用于非周期离散信号% 生成非周期离散信号 n 0:100; x (0.9).^n; % 指数衰减序列2. MATLAB实现与可视化对比2.1 周期连续信号的FS分析让我们从一个简单的周期信号开始——方波信号% 方波信号的FS分析 T 2; % 周期(秒) f0 1/T; % 基频(Hz) t linspace(0, 3*T, 1000); % 时间向量 square_wave square(2*pi*f0*t); % 生成方波 % 计算FS系数(前20个谐波) N 20; k 1:N; ak (2./(pi*k)) .* sin(pi*k/2) .* (mod(k,2)1); % FS系数 bk zeros(size(k)); % 方波只有奇次谐波 % 重构信号 reconstructed zeros(size(t)); for n 1:length(k) reconstructed reconstructed ak(n)*cos(2*pi*n*f0*t) bk(n)*sin(2*pi*n*f0*t); end % 绘制结果 figure; subplot(2,1,1); plot(t, square_wave); title(原始方波信号); subplot(2,1,2); plot(t, reconstructed); title(FS重构信号);关键观察点FS系数ak表示各频率分量的幅度方波只包含奇次谐波随着谐波数量增加重构信号越来越接近原始信号2.2 非周期连续信号的FT分析考虑一个非周期信号——高斯脉冲的FT分析% 高斯脉冲的FT分析 t linspace(-5, 5, 1000); sigma 0.5; gaussian exp(-t.^2/(2*sigma^2)); % 高斯脉冲 % 计算FT(使用FFT近似) Fs 1000/(t(end)-t(1)); % 采样频率 N length(gaussian); f (-N/2:N/2-1)*(Fs/N); % 频率向量 ft_gaussian fftshift(fft(gaussian))/N; % FFT并归一化 % 理论FT结果 theory_ft sigma*sqrt(2*pi)*exp(-(2*pi*sigma*f).^2/2); % 绘制结果 figure; subplot(2,1,1); plot(t, gaussian); title(高斯脉冲); subplot(2,1,2); plot(f, abs(ft_gaussian), f, theory_ft, --); legend(FFT结果, 理论FT); title(幅度谱);重要发现高斯脉冲的FT仍然是高斯函数FFT可以很好地近似连续FT频域宽度与时域宽度成反比(不确定性原理)2.3 周期离散信号的DFS分析分析一个周期离散信号——周期方波序列% 周期离散信号的DFS分析 N 20; % 周期 n 0:3*N-1; % 3个周期 x mod(n, N) N/2; % 周期方波序列 % 计算DFS系数 k 0:N-1; Xk zeros(size(k)); for kk k Xk(kk1) sum(x(1:N).*exp(-1j*2*pi*kk*(0:N-1)/N)); end % 重构信号 x_recon zeros(size(n)); for nn n x_recon(nn1) (1/N)*sum(Xk.*exp(1j*2*pi*(0:N-1)*nn/N)); end % 绘制结果 figure; subplot(3,1,1); stem(n, x); title(原始周期序列); subplot(3,1,2); stem(k, abs(Xk)); title(DFS幅度谱); subplot(3,1,3); stem(n, real(x_recon)); title(DFS重构信号);关键点DFS频谱也是离散且周期的重构信号与原始信号完全一致(无信息丢失)离散信号的频谱是周期的(2π周期)2.4 非周期离散信号的DTFT分析分析一个非周期离散信号——指数衰减序列% 非周期离散信号的DTFT分析 n 0:100; a 0.9; x a.^n; % 指数衰减序列 % 计算DTFT w linspace(-pi, pi, 1000); % 数字频率 X zeros(size(w)); for k 1:length(w) X(k) sum(x .* exp(-1j*w(k)*n)); end % 理论DTFT结果 theory_X 1./(1 - a*exp(-1j*w)); % 绘制结果 figure; subplot(2,1,1); stem(n, x); title(指数衰减序列); subplot(2,1,2); plot(w, abs(X), w, abs(theory_X), --); legend(数值计算, 理论结果); title(DTFT幅度谱);重要特性DTFT频谱是连续的且2π周期的对于这个简单序列DTFT有解析解数值计算结果与理论解完全吻合3. 工程应用场景与选择指南不同的傅里叶分析方法适用于不同的工程场景。下面通过具体案例说明如何选择适当的分析方法。3.1 音频信号处理音频信号通常是连续时间信号但在数字系统中我们处理的是其采样后的离散版本。典型处理流程模拟信号通过抗混叠滤波器ADC采样得到离散信号使用DFT/FFT进行频谱分析数字信号处理DAC重建模拟信号% 音频信号频谱分析示例 [audio, Fs] audioread(sample.wav); % 读取音频文件 N length(audio); f (0:N-1)*(Fs/N); % 频率轴 % 计算频谱 spectrum abs(fft(audio)); % 绘制频谱 figure; plot(f(1:N/2), spectrum(1:N/2)); xlabel(频率(Hz)); ylabel(幅度); title(音频信号频谱);注意实际音频处理中通常使用短时傅里叶变换(STFT)来分析时变频谱特性。3.2 通信系统中的调制分析在通信系统中FS和FT常用于分析调制信号的频谱特性。AM调制信号分析% AM调制信号分析 fc 1000; % 载波频率(Hz) fm 100; % 调制频率(Hz) Fs 10000; % 采样率 t 0:1/Fs:1; % 1秒时间 % 生成AM信号 m 0.5; % 调制深度 carrier cos(2*pi*fc*t); message cos(2*pi*fm*t); AM (1 m*message).*carrier; % 计算频谱 N length(AM); f (-N/2:N/2-1)*(Fs/N); AM_spectrum abs(fftshift(fft(AM)))/N; % 绘制结果 figure; subplot(2,1,1); plot(t, AM); title(AM调制信号); subplot(2,1,2); plot(f, AM_spectrum); xlim([-1500 1500]); title(AM信号频谱);频谱特征载波频率fc处的峰值边带(fc±fm)处的峰值频谱对称性3.3 数字滤波器设计DTFT在数字滤波器设计和分析中起着核心作用。FIR滤波器频率响应分析% FIR滤波器设计与频率响应分析 N 50; % 滤波器阶数 fc 0.2; % 截止频率(归一化) b fir1(N, fc); % 设计FIR滤波器 % 计算频率响应 [w, H] freqz(b, 1, 1024); % 绘制结果 figure; subplot(2,1,1); stem(b); title(滤波器系数); subplot(2,1,2); plot(w/pi, 20*log10(abs(H))); xlabel(归一化频率(\pi rad/sample)); ylabel(幅度(dB)); title(频率响应);关键概念频率响应就是滤波器冲激响应的DTFT归一化频率(1对应Nyquist频率)幅度响应和相位响应4. 常见问题与实用技巧在实际应用中工程师经常会遇到一些典型问题和挑战。本节提供一些实用解决方案。4.1 频谱泄漏与窗函数当分析有限长度信号时频谱泄漏是一个常见问题。使用适当的窗函数可以减少泄漏。% 频谱泄漏与窗函数比较 N 100; n 0:N-1; f0 0.3; % 非整数倍周期 x cos(2*pi*f0*n); % 不加窗和加汉宁窗的频谱比较 X_nowin abs(fft(x)); win hann(N); X_win abs(fft(x.*win)); % 绘制结果 figure; subplot(2,1,1); stem(n, x); title(原始信号); subplot(2,1,2); plot(0:N-1, 20*log10(X_nowin), b, ... 0:N-1, 20*log10(X_win), r); legend(无窗, 汉宁窗); title(频谱对比);窗函数选择指南窗函数主瓣宽度旁瓣衰减适用场景矩形窗窄差(-13dB)瞬态信号精确频率测量汉宁窗中等好(-31dB)一般频谱分析汉明窗中等较好(-42dB)音频处理布莱克曼窗宽很好(-58dB)要求高动态范围表2常用窗函数特性比较4.2 频率分辨率与补零频率分辨率取决于信号实际长度补零可以提高频谱显示的平滑度但不增加实际分辨率。% 频率分辨率与补零演示 N 32; % 实际数据长度 n 0:N-1; x cos(2*pi*0.2*n) 0.5*cos(2*pi*0.3*n); % 不同补零长度的FFT X1 abs(fft(x)); X2 abs(fft(x, 64)); X3 abs(fft(x, 128)); % 绘制结果 figure; subplot(3,1,1); stem(0:N-1, X1); title(无补零); subplot(3,1,2); stem(0:63, X2); title(补零到64点); subplot(3,1,3); stem(0:127, X3); title(补零到128点);重要结论频率分辨率ΔfFs/N补零不能提高Δf补零使频谱显示更平滑4.3 实时频谱分析实现在许多应用中需要实时计算和显示信号频谱。下面是一个简单的实时频谱分析框架。% 简易实时频谱分析框架 Fs 44100; % 采样率 frameSize 1024; % 帧大小 spectrumLength frameSize/2; % 初始化图形 figure; h plot(zeros(1, spectrumLength)); ylim([0 100]); xlabel(频率(bin)); ylabel(幅度); title(实时频谱分析); % 模拟实时处理 while ishandle(h) % 模拟获取一帧音频数据(实际应用中替换为真实音频输入) frame randn(1, frameSize); % 使用噪声作为示例 % 加窗并计算FFT win hann(frameSize); frame_win frame .* win; spectrum abs(fft(frame_win)); % 更新显示 set(h, YData, 20*log10(spectrum(1:spectrumLength))); drawnow; end优化建议使用重叠-保留法减少帧间突变采用移动平均平滑频谱显示对于嵌入式系统优化FFT计算效率5. 高级应用与扩展5.1 时频分析短时傅里叶变换对于非平稳信号单纯的傅里叶分析不足以描述其频率特性随时间的变化。短时傅里叶变换(STFT)提供了时频联合分析的能力。% STFT示例啁啾信号分析 Fs 1000; t 0:1/Fs:1; f0 50; f1 200; x chirp(t, f0, 1, f1); % 线性啁啾信号 % 计算STFT window hann(128); noverlap 120; nfft 256; [S, f, t_stft] spectrogram(x, window, noverlap, nfft, Fs); % 绘制时频谱 figure; imagesc(t_stft, f, 20*log10(abs(S))); axis xy; colorbar; xlabel(时间(s)); ylabel(频率(Hz)); title(啁啾信号的STFT);STFT参数选择窗长度决定时间/频率分辨率折衷重叠率通常50-75%减少帧间突变FFT点数通常≥窗长度可补零5.2 多速率信号处理在多速率系统中不同采样率转换时需要特别注意频域特性的变化。% 采样率转换与频谱变化 Fs_orig 1000; % 原始采样率 t_orig 0:1/Fs_orig:1; x_orig sin(2*pi*100*t_orig) 0.5*sin(2*pi*300*t_orig); % 降采样 M 2; % 降采样因子 Fs_new Fs_orig/M; x_down x_orig(1:M:end); % 绘制频谱比较 figure; subplot(2,1,1); plot((0:length(x_orig)-1)/length(x_orig)*Fs_orig, abs(fft(x_orig))); title(原始信号频谱); xlim([0 Fs_orig/2]); subplot(2,1,2); plot((0:length(x_down)-1)/length(x_down)*Fs_new, abs(fft(x_down))); title(降采样后频谱); xlim([0 Fs_new/2]);抗混叠处理降采样前必须进行低通滤波截止频率≤新Nyquist频率MATLAB中可用decimate函数自动处理5.3 倒谱分析倒谱分析在语音处理和故障诊断中有广泛应用可以分离激励源和系统响应。% 倒谱分析示例 Fs 8000; [y, Fs] audioread(speech.wav); % 读取语音样本 % 计算实倒谱 cepstrum ifft(log(abs(fft(y)))); % 绘制结果 figure; subplot(2,1,1); plot(y); title(语音信号); subplot(2,1,2); plot((0:length(cepstrum)-1)/Fs, abs(cepstrum)); xlabel(倒频率(s)); title(倒谱); xlim([0 0.02]); % 显示低倒频率部分倒谱应用基音周期检测声道响应分析机械故障特征提取