使用汉克矩阵构造信号，并通过差分谱法进行奇异值分解(SVD)去噪。该算法特别适用于处理含噪的周期信号或准周期信号。
function svd_denoising_hankel()
% 清空工作区
clc;
close all;
clear;
%% 1. 生成模拟信号
fs = 1000; % 采样频率 (Hz)
t = 0:1/fs:1-1/fs; % 时间向量 (1秒)
N = length(t); % 信号长度
% 原始信号: 多个正弦波叠加
f1 = 50; % 基频 (Hz)
f2 = 120; % 二次谐波
f3 = 80; % 另一个频率成分
signal = 1.5*sin(2*pi*f1*t) + 0.8*cos(2*pi*f2*t) + 1.2*sin(2*pi*f3*t);
% 添加高斯白噪声
SNR = 5; % 信噪比 (dB)
noisy_signal = awgn(signal, SNR, 'measured');
%% 2. 构造汉克矩阵
L = floor(N/2); % 汉克矩阵行数 (经验值)
K = N - L + 1; % 汉克矩阵列数
% 构建汉克矩阵
H = zeros(L, K);
for i = 1:L
H(i, :) = noisy_signal(i:i+K-1);
end
%% 3. 奇异值分解(SVD)
[U, S, V] = svd(H, 'econ');
singular_values = diag(S); % 提取奇异值
%% 4. 差分谱分析确定有效秩
% 计算差分谱
diff_spectrum = abs(diff(singular_values));
% 找到差分谱的峰值点
[pks, locs] = findpeaks(diff_spectrum, 'MinPeakHeight', 0.1*max(diff_spectrum));
% 确定有效秩 (保留前r个奇异值)
if isempty(locs)
r = 1; % 如果没有找到峰值，至少保留一个奇异值
else
[~, idx] = max(pks); % 找到最大峰值
r = locs(idx); % 有效秩
end
% 确保有效秩合理
r = min(r, length(singular_values)-1);
r = max(r, 1);
%% 5. 信号重构
% 保留前r个奇异值
S_denoised = S;
S_denoised(r+1:end, r+1:end) = 0;
% 重构汉克矩阵
H_denoised = U * S_denoised * V';
% 从汉克矩阵重构信号 (反对角线平均)
denoised_signal = zeros(1, N);
count = zeros(1, N);
for i = 1:L
for j = 1:K
idx = i + j - 1; % 反对角线索引
denoised_signal(idx) = denoised_signal(idx) + H_denoised(i, j);
count(idx) = count(idx) + 1;
end
end
% 计算平均值
denoised_signal = denoised_signal ./ count;
%% 6. 性能评估
% 计算信噪比改进
noise_before = noisy_signal - signal;
noise_after = denoised_signal - signal;
SNR_before = 10*log10(var(signal)/var(noise_before));
SNR_after = 10*log10(var(signal)/var(noise_after));
improvement = SNR_after - SNR_before;
% 计算均方根误差
RMSE_before = sqrt(mean(noise_before.^2));
RMSE_after = sqrt(mean(noise_after.^2));
%% 7. 结果可视化
% 创建主图窗
fig = figure('Name', '基于汉克矩阵和差分谱法的SVD信号去噪', ...
'Position', [100, 100, 1200, 800], ...
'Color', [0.95, 0.95, 0.95]);
% 信号对比图
subplot(3, 2, [1, 2]);
plot(t, signal, 'b', 'LineWidth', 1.8, 'DisplayName', '原始信号');
hold on;
plot(t, noisy_signal, 'Color', [0.7, 0.7, 0.7], 'LineWidth', 1, 'DisplayName', '含噪信号');
plot(t, denoised_signal, 'r', 'LineWidth', 1.5, 'DisplayName', '去噪信号');
hold off;
title(sprintf('信号对比 (SNR提升: %.2f dB)', improvement), 'FontSize', 12);
xlabel('时间 (s)');
ylabel('幅值');
legend('Location', 'best');
grid on;
xlim([0, 0.2]); % 显示前200ms
% 频谱分析
subplot(3, 2, 3);
[freq, P1] = calc_spectrum(signal, fs);
[~, P2] = calc_spectrum(noisy_signal, fs);
[~, P3] = calc_spectrum(denoised_signal, fs);
plot(freq, P1, 'b', 'LineWidth', 1.5);
hold on;
plot(freq, P2, 'Color', [0.7, 0.7, 0.7], 'LineWidth', 1);
plot(freq, P3, 'r', 'LineWidth', 1.2);
hold off;
title('频谱分析');
xlabel('频率 (Hz)');
ylabel('功率谱密度');
grid on;
xlim([0, 200]);
% 奇异值分布
subplot(3, 2, 4);
semilogy(singular_values, 'bo-', 'LineWidth', 1.5, 'MarkerSize', 4);
hold on;
semilogy(1:r, singular_values(1:r), 'ro', 'MarkerSize', 6, 'LineWidth', 1.5);
line([r+0.5, r+0.5], ylim, 'Color', 'k', 'LineStyle', '--', 'LineWidth', 1.2);
hold off;
title(sprintf('奇异值分布 (有效秩: r = %d)', r));
xlabel('奇异值索引');
ylabel('奇异值 (对数尺度)');
grid on;
% 差分谱分析
subplot(3, 2, 5);
plot(diff_spectrum, 'go-', 'LineWidth', 1.5, 'MarkerSize', 4);
hold on;
plot(locs, pks, 'r*', 'MarkerSize', 10, 'LineWidth', 1.8);
plot([r, r], [0, max(diff_spectrum)], 'k--', 'LineWidth', 1.2);
hold off;
title('奇异值差分谱');
xlabel('索引');
ylabel('差分值');
grid on;
xlim([1, length(diff_spectrum)]);
% 残差分析
subplot(3, 2, 6);
plot(t, noise_before, 'Color', [0.7, 0.7, 0.7], 'LineWidth', 1);
hold on;
plot(t, noise_after, 'r', 'LineWidth', 1.2);
hold off;
title(sprintf('噪声残差 (RMSE: %.4f -> %.4f)', RMSE_before, RMSE_after));
xlabel('时间 (s)');
ylabel('残差');
legend('去噪前残差', '去噪后残差');
grid on;
xlim([0, 0.2]); % 显示前200ms
% 添加信息面板
info_str = sprintf(['算法参数:\n' ...
' 信号长度: %d\n' ...
' 汉克矩阵: %d×%d\n' ...
' 信噪比(输入): %.2f dB\n' ...
' 信噪比(输出): %.2f dB\n' ...
' 信噪比提升: %.2f dB\n' ...
' RMSE减少: %.2f%%'], ...
N, L, K, SNR_before, SNR_after, improvement, ...
100*(RMSE_before-RMSE_after)/RMSE_before);
annotation('textbox', [0.1, 0.01, 0.8, 0.08], ...
'String', info_str, ...
'FitBoxToText', 'on', ...
'BackgroundColor', [0.98, 0.98, 0.98], ...
'EdgeColor', [0.8, 0.8, 0.8], ...
'FontSize', 10);
%% 辅助函数：计算频谱
function [f, P] = calc_spectrum(x, fs)
n = length(x);
y = fft(x);
P2 = abs(y/n);
P1 = P2(1:floor(n/2)+1);
P1(2:end-1) = 2*P1(2:end-1);
f = fs*(0:(n/2))/n;
f = f(:);
P = P1(:);
end
end
算法原理与实现步骤
1. 汉克矩阵构造
汉克矩阵是一种特殊矩阵，其反对角线上的元素相同：
H = zeros(L, K);
for i = 1:L
H(i, :) = noisy_signal(i:i+K-1);
end
其中：
L是矩阵行数（通常取信号长度的一半）K = N - L + 1是矩阵列数N是信号长度
汉克矩阵的构造将一维信号转换为二维矩阵，为后续的SVD分解做准备。
2. 奇异值分解(SVD)
对汉克矩阵进行SVD分解：
[U, S, V] = svd(H, 'econ');
singular_values = diag(S);
SVD将矩阵分解为：
H = U * S * V'
其中：
U是左奇异向量矩阵S是对角奇异值矩阵V是右奇异向量矩阵
3. 差分谱分析
差分谱是确定有效秩的关键：
diff_spectrum = abs(diff(singular_values));
[pks, locs] = findpeaks(diff_spectrum, 'MinPeakHeight', 0.1*max(diff_spectrum));
奇异值差分谱定义为：
Δσ(i) = |σ(i) - σ(i+1)|
其中σ(i)是第i个奇异值。
有效秩r由差分谱的最大峰值位置确定：
[~, idx] = max(pks);
r = locs(idx);
4. 信号重构
基于有效秩重构信号：
S_denoised = S;
S_denoised(r+1:end, r+1:end) = 0;
H_denoised = U * S_denoised * V';
从重构的汉克矩阵恢复一维信号（反对角线平均法）：
denoised_signal = zeros(1, N);
count = zeros(1, N);
for i = 1:L
for j = 1:K
idx = i + j - 1;
denoised_signal(idx) = denoised_signal(idx) + H_denoised(i, j);
count(idx) = count(idx) + 1;
end
end
denoised_signal = denoised_signal ./ count;
参考代码 采用汉克矩阵构造信号，使用差分谱法进行SVD去噪 www.youwenfan.com/contentalf/96370.html
扩展与改进
自适应汉克矩阵大小：
% 基于信号自相关确定L [acf, ~] = xcorr(noisy_signal, 'coeff'); acf = acf(length(noisy_signal):end); [~, lags] = findpeaks(acf, 'MinPeakHeight', 0.5); if ~isempty(lags) L = min(lags(1), floor(N/2)); else L = floor(N/2); end多通道信号处理：
% 对于多通道信号，构建块汉克矩阵 H_multi = []; for ch = 1:nChannels H_ch = hankel_matrix(signal_ch(:, ch), L); H_multi = [H_multi; H_ch]; end实时处理实现：
% 滑动窗口处理 window_size = 1000; step_size = 100; for start_idx = 1:step_size:(N-window_size) segment = noisy_signal(start_idx:start_idx+window_size-1); denoised_segment = svd_denoise(segment); % 拼接处理结果 end结合小波阈值：
% SVD去噪后的小波阈值处理 denoised_signal_wavelet = wden(denoised_signal, 'modwtsqtwolog', 's', 'mln', 5, 'db4');
这个实现提供了基于汉克矩阵和差分谱法的完整SVD去噪流程，通过可视化展示各步骤结果，便于理解和分析算法性能。
