Matlab STFT 频谱图和幅度图优化

使用 Matlab 进行 STFT 分析，优化频谱图和幅度图

本文将使用 Matlab 代码演示如何进行短时傅里叶变换 (STFT) 分析，并提供优化频谱图精度和将幅度图转换为折线图的技巧。

代码示例

ns = 1024;
t = linspace(0, 1024, ns);
f = 50;
w = 2 * pi * f;

Voltage_amplitude = 10000; % 10kV
sinewave_voltage = Voltage_amplitude * sin(w * t);

Current_amplitude = 120; % 120A
sinewave_current = Current_amplitude * sin(w * t);

alpha1_h = 1;
alpha3_h = 0.015;
alpha5_h = 0.01;
alpha7_h = 0.09;

harmonics_current = alpha1_h * sinewave_current + alpha3_h * sin(3 * w * t) + alpha5_h * sin(5 * w * t) + alpha7_h * sin(7 * w * t);

signal = harmonics_current;

win = 256; % 窗口大小
hop = 64; % 跳跃大小
nfft = 1024; % FFT点数
fs = 1000; % 采样频率，单位是赫兹

% 计算STFT
[STFT, f, t] = stft(harmonics_current, win, hop, nfft, fs);

% 计算幅度
M = abs(STFT);

% 计算功率谱密度（用于绘制频谱图）
S = 20 * log10(M);

% 绘制频谱图
figure;
imagesc(t, f, S);
axis xy;
xlabel('时间（秒）');
ylabel('频率（赫兹）');
title('频谱图');
colorbar;

% 绘制幅度图
figure;
plot(t, M);
xlabel('时间（秒）');
ylabel('幅度');
title('幅度图');

function [S, f, t] = stft(x, win, hop, nfft, fs)
% 计算窗口长度
win_length = length(win);
% 计算信号长度
signal_length = length(x);
% 计算帧数
num_frames = floor((signal_length - win_length) / hop) + 1;
% 初始化STFT矩阵
S = zeros(nfft, num_frames);
% 对每一帧进行FFT
for k = 1:num_frames
% 提取当前帧
frame = x(((k - 1) * hop + 1):((k - 1) * hop + win_length));
% 应用窗口函数
windowed_frame = frame .* win;
% 计算FFT
S(:, k) = fft(windowed_frame, nfft);
end
% 计算频率向量
f = (0:(nfft / 2 - 1)) * fs / nfft;
% 计算时间向量
t = (0:(num_frames - 1)) * hop / fs;
% 只保留正频率部分
S = S(1:nfft / 2, :);
end

优化频谱图精度

• 增加 nfft 的值以增加 FFT 点数，从而提高频谱分辨率。
• 降低 hop 的值以增加帧数，从而提高时间分辨率。

例如，在代码中，我们将 nfft 调整为 1024，hop 调整为 64。

使用折线图绘制幅度图

• 可以使用 plot 函数绘制频率随时间变化的曲线，而不是使用 imagesc 函数。

代码示例中，我们使用以下代码绘制幅度图：

figure;
plot(t, M);
xlabel('时间（秒）');
ylabel('幅度');
title('幅度图');

总结

通过以上优化，我们可以获得更精确的频谱图，并将幅度图以折线图的形式展示出来，更直观地展现信号的时频特性。