快速傅里叶变换(FFT)的应用
1. 信号处理基础概念
%% FFT基础与应用示例
clear; close all; clc;
% 生成测试信号
fs = 1000; % 采样频率 1000Hz
t = 0:1/fs:1-1/fs; % 时间向量 (1秒)
N = length(t); % 信号长度
% 创建包含多个频率成分的信号
f1 = 50; % 50Hz
f2 = 120; % 120Hz
f3 = 300; % 300Hz
signal = 0.7*sin(2*pi*f1*t) + sin(2*pi*f2*t) + 0.3*sin(2*pi*f3*t);
noise = 0.5*randn(size(t)); % 高斯白噪声
x = signal + noise; % 含噪声的信号
% 绘制时域信号
figure('Position', [100, 100, 1200, 800]);
subplot(2,2,1);
plot(t, x);
title('时域信号');
xlabel('时间 (s)');
ylabel('幅度');
grid on;
xlim([0 0.1]); % 只看前0.1秒
subplot(2,2,2);
plot(t, signal);
title('纯净信号');
xlabel('时间 (s)');
ylabel('幅度');
grid on;
xlim([0 0.1]);
2. FFT频谱分析完整实现
%% FFT频谱分析函数
function [f, P1, P2] = computeSpectrum(x, fs, window_type)
% 计算信号的频谱
% 输入:
% x - 输入信号
% fs - 采样频率
% window_type - 窗函数类型 ('hann', 'hamming', 'blackman', 'rectangular')
% 输出:
% f - 频率向量
% P1 - 单边幅度谱
% P2 - 双边幅度谱
N = length(x);
% 应用窗函数
switch window_type
case 'hann'
w = hann(N)';
case 'hamming'
w = hamming(N)';
case 'blackman'
w = blackman(N)';
otherwise
w = ones(1, N); % 矩形窗
end
x_windowed = x .* w;
% 计算FFT
X = fft(x_windowed);
% 计算双边谱
P2 = abs(X/N);
% 计算单边谱
P1 = P2(1:floor(N/2)+1);
P1(2:end-1) = 2*P1(2:end-1);
% 频率向量
f = fs*(0:(N/2))/N;
end
%% 频谱分析应用
% 使用不同窗函数计算频谱
[f_rect, P1_rect, P2_rect] = computeSpectrum(x, fs, 'rectangular');
[f_hann, P1_hann, P2_hann] = computeSpectrum(x, fs, 'hann');
[f_ham, P1_ham, P2_ham] = computeSpectrum(x, fs, 'hamming');
% 绘制频谱比较
figure('Position', [100, 100, 1400, 900]);
% 单边幅度谱比较
subplot(2,3,1);
plot(f_rect, P1_rect, 'b-', 'LineWidth', 1.5);
title('矩形窗 - 单边幅度谱');
xlabel('频率 (Hz)');
ylabel('|P1(f)|');
grid on;
xlim([0 500]);
subplot(2,3,2);
plot(f_hann, P1_hann, 'r-', 'LineWidth', 1.5);
title('汉宁窗 - 单边幅度谱');
xlabel('频率 (Hz)');
ylabel('|P1(f)|');
grid on;
xlim([0 500]);
subplot(2,3,3);
plot(f_ham, P1_ham, 'g-', 'LineWidth', 1.5);
title('汉明窗 - 单边幅度谱');
xlabel('频率 (Hz)');
ylabel('|P1(f)|');
grid on;
xlim([0 500]);
% 双边幅度谱比较
subplot(2,3,4);
f_full = fs*(0:N-1)/N;
plot(f_full, P2_rect, 'b-', 'LineWidth', 1.5);
title('矩形窗 - 双边幅度谱');
xlabel('频率 (Hz)');
ylabel('|P2(f)|');
grid on;
xlim([0 fs]);
subplot(2,3,5);
plot(f_full, P2_hann, 'r-', 'LineWidth', 1.5);
title('汉宁窗 - 双边幅度谱');
xlabel('频率 (Hz)');
ylabel('|P2(f)|');
grid on;
xlim([0 fs]);
subplot(2,3,6);
plot(f_full, P2_ham, 'g-', 'LineWidth', 1.5);
title('汉明窗 - 双边幅度谱');
xlabel('频率 (Hz)');
ylabel('|P2(f)|');
grid on;
xlim([0 fs]);
3. 功率谱密度估计
%% 功率谱密度分析
function [f, psd] = computePSD(x, fs, method)
% 计算功率谱密度
% method: 'periodogram', 'welch'
N = length(x);
switch method
case 'periodogram'
% 周期图法
[psd, f] = periodogram(x, hann(N), N, fs);
case 'welch'
% Welch方法
window = hann(256);
noverlap = 128;
nfft = 512;
[psd, f] = pwelch(x, window, noverlap, nfft, fs);
end
end
% 计算不同方法的PSD
[f_per, psd_per] = computePSD(x, fs, 'periodogram');
[f_welch, psd_welch] = computePSD(x, fs, 'welch');
% 绘制功率谱密度
figure('Position', [100, 100, 1200, 500]);
subplot(1,2,1);
plot(f_per, 10*log10(psd_per), 'b-', 'LineWidth', 1.5);
title('周期图法 - 功率谱密度');
xlabel('频率 (Hz)');
ylabel('功率谱密度 (dB/Hz)');
grid on;
xlim([0 500]);
subplot(1,2,2);
plot(f_welch, 10*log10(psd_welch), 'r-', 'LineWidth', 1.5);
title('Welch方法 - 功率谱密度');
xlabel('频率 (Hz)');
ylabel('功率谱密度 (dB/Hz)');
grid on;
xlim([0 500]);
4. 实际应用案例
%% 应用案例1: 音频信号分析
% 生成模拟音频信号（钢琴音符）
fs_audio = 44100; % 音频采样率
t_audio = 0:1/fs_audio:1-1/fs_audio;
% C大调和弦: C4(261.63Hz), E4(329.63Hz), G4(392.00Hz)
f_C4 = 261.63;
f_E4 = 329.63;
f_G4 = 392.00;
chord = sin(2*pi*f_C4*t_audio) + 0.8*sin(2*pi*f_E4*t_audio) + 0.6*sin(2*pi*f_G4*t_audio);
chord = chord/max(abs(chord)); % 归一化
% 频谱分析
[f_chord, P1_chord] = computeSpectrum(chord, fs_audio, 'hann');
figure('Position', [100, 100, 1000, 600]);
subplot(2,1,1);
plot(t_audio(1:2000), chord(1:2000));
title('C大调和弦时域波形');
xlabel('时间 (s)');
ylabel('幅度');
grid on;
subplot(2,1,2);
plot(f_chord, 20*log10(P1_chord), 'b-', 'LineWidth', 1.5);
title('和弦频谱');
xlabel('频率 (Hz)');
ylabel('幅度 (dB)');
grid on;
xlim([200 500]);
% 标记基频
hold on;
plot([f_C4, f_E4, f_G4], [0, 0, 0], 'ro', 'MarkerSize', 8, 'LineWidth', 2);
legend('频谱', '基频位置', 'Location', 'best');
5. 信号调制与解调分析
%% 应用案例2: 调制信号分析
% AM调制信号分析
fc = 100; % 载波频率
fm = 10; % 调制频率
A = 1; % 载波幅度
m = 0.8; % 调制深度
t_mod = 0:1/fs:1-1/fs;
carrier = A * sin(2*pi*fc*t_mod);
modulating = sin(2*pi*fm*t_mod);
am_signal = A * (1 + m*modulating) .* sin(2*pi*fc*t_mod);
% 频谱分析
[f_am, P1_am] = computeSpectrum(am_signal, fs, 'hann');
figure('Position', [100, 100, 1200, 800]);
% 时域信号
subplot(3,1,1);
plot(t_mod, modulating, 'g-', 'LineWidth', 1.5);
hold on;
plot(t_mod, am_signal, 'b-', 'LineWidth', 1);
title('AM调制信号');
xlabel('时间 (s)');
ylabel('幅度');
legend('调制信号', '已调信号', 'Location', 'best');
grid on;
xlim([0 0.5]);
% 载波信号
subplot(3,1,2);
plot(t_mod(1:500), carrier(1:500), 'r-', 'LineWidth', 1.5);
title('载波信号');
xlabel('时间 (s)');
ylabel('幅度');
grid on;
% 频谱
subplot(3,1,3);
plot(f_am, 20*log10(P1_am), 'b-', 'LineWidth', 1.5);
title('AM信号频谱');
xlabel('频率 (Hz)');
ylabel('幅度 (dB)');
grid on;
xlim([0 200]);
% 标记载波和边带频率
hold on;
plot([fc-fm, fc, fc+fm], [-20, 0, -20], 'ro', 'MarkerSize', 8, 'LineWidth', 2);
legend('频谱', '载波和边带', 'Location', 'best');
6. 频谱特征提取
%% 频谱特征提取函数
function features = extractSpectralFeatures(x, fs)
% 提取频谱特征
[f, P1] = computeSpectrum(x, fs, 'hann');
features = struct();
% 找到主要频率峰值
[peaks, locs] = findpeaks(P1, 'SortStr', 'descend', 'NPeaks', 5);
features.peak_frequencies = f(locs(1:min(3, length(locs))));
features.peak_magnitudes = peaks(1:min(3, length(peaks)));
% 频谱质心
features.spectral_centroid = sum(f .* P1) / sum(P1);
% 频谱带宽
features.spectral_bandwidth = sqrt(sum((f - features.spectral_centroid).^2 .* P1) / sum(P1));
% 频谱滚降点 (85%)
total_power = sum(P1);
cumulative_power = cumsum(P1);
rolloff_idx = find(cumulative_power >= 0.85 * total_power, 1);
features.spectral_rolloff = f(rolloff_idx);
% 频谱平坦度
features.spectral_flatness = exp(mean(log(P1 + eps))) / (mean(P1) + eps);
% 频谱熵
P_norm = P1 / sum(P1);
features.spectral_entropy = -sum(P_norm .* log2(P_norm + eps));
end
% 特征提取演示
features_signal = extractSpectralFeatures(signal, fs);
features_noisy = extractSpectralFeatures(x, fs);
fprintf('=== 频谱特征比较 ===\n');
fprintf('纯净信号特征:\n');
fprintf(' 频谱质心: %.2f Hz\n', features_signal.spectral_centroid);
fprintf(' 频谱带宽: %.2f Hz\n', features_signal.spectral_bandwidth);
fprintf(' 频谱滚降: %.2f Hz\n', features_signal.spectral_rolloff);
fprintf(' 频谱平坦度: %.4f\n', features_signal.spectral_flatness);
fprintf(' 频谱熵: %.4f\n', features_signal.spectral_entropy);
fprintf('\n含噪声信号特征:\n');
fprintf(' 频谱质心: %.2f Hz\n', features_noisy.spectral_centroid);
fprintf(' 频谱带宽: %.2f Hz\n', features_noisy.spectral_bandwidth);
fprintf(' 频谱滚降: %.2f Hz\n', features_noisy.spectral_rolloff);
fprintf(' 频谱平坦度: %.4f\n', features_noisy.spectral_flatness);
fprintf(' 频谱熵: %.4f\n', features_noisy.spectral_entropy);
7. 实时频谱分析模拟
%% 实时频谱分析模拟
function realTimeSpectrumDemo()
% 模拟实时频谱分析
fs = 1000;
t_chunk = 0.1; % 每段分析0.1秒
samples_per_chunk = round(fs * t_chunk);
figure('Position', [200, 200, 1000, 600]);
for i = 1:20
% 生成随时间变化的信号
t = (0:samples_per_chunk-1)/fs + (i-1)*t_chunk;
% 频率随时间变化
f_current = 50 + 10*sin(2*pi*0.5*t);
signal_chunk = sin(2*pi*f_current.*t) + 0.3*randn(size(t));
% 计算频谱
[f, P1] = computeSpectrum(signal_chunk, fs, 'hann');
% 更新显示
subplot(2,1,1);
plot(t, signal_chunk, 'b-', 'LineWidth', 1);
title(sprintf('时域信号 (块 %d)', i));
xlabel('时间 (s)');
ylabel('幅度');
grid on;
subplot(2,1,2);
plot(f, 20*log10(P1), 'r-', 'LineWidth', 1.5);
title('频谱');
xlabel('频率 (Hz)');
ylabel('幅度 (dB)');
grid on;
xlim([0 200]);
ylim([-60 0]);
drawnow;
pause(0.5); % 模拟实时延迟
end
end
% 运行实时频谱分析演示
% realTimeSpectrumDemo();
参考代码 快速傅里叶变换的应用以及傅里叶变换求频谱 www.youwenfan.com/contentcnk/79398.html
知识点总结
FFT应用领域：
- 音频处理 - 音乐分析、语音识别
- 通信系统 - 调制解调、频谱监测
- 振动分析 - 机械故障诊断
- 图像处理 - 频域滤波
- 生物医学 - EEG、ECG信号分析
频谱分析要点：
- 采样定理 - 采样频率 > 2倍最高频率
- 窗函数选择 - 减少频谱泄漏
- 频率分辨率 - Δf = fs/N
- 幅度校正 - 单边谱需要乘以2（除直流分量）
性能优化技巧：
- 使用适当的窗函数减少频谱泄漏
- 零填充提高频率显示分辨率
- 平均处理降低噪声影响
- 选择合适的FFT长度（2的幂次）
浙公网安备 33010602011771号