CLEAN 算法仿真：雷达信号杂波抑制

CLEAN 算法仿真程序，用于雷达信号中的杂波抑制（如旁瓣或强散射点干扰）。

CLEAN 算法本质是一种迭代反卷积技术，通过识别和减去最强散射分量来'清理'信号。

CLEAN 算法核心原理

CLEAN 算法的核心思想非常直观，其处理流程如下：

1. 输入带噪信号
2. 在测量谱中寻找最大峰值
3. 判断峰值强度是否大于阈值
4. 是：生成并存储该分量（位置、幅度），从原始信号中减去该分量模型，更新残差信号
5. 否：结束循环
6. 重构'干净'信号（所有分量合成）
7. 输出高分辨率的干净信号

MATLAB 仿真程序

%% CLEAN 算法用于杂波抑制 - 完整仿真程序
clear; clc; close all;

%% 1. 参数设置
N = 256; % 信号长度
f0 = 0.2; % 主信号频率（归一化）
f_clutter = [0.1, 0.25, 0.35]; % 杂波频率
A_clutter = [0.5, 0.8, 1.2]; % 杂波幅度（强于目标）
SNR_dB = 15; % 信噪比 (dB)
K = 20; % CLEAN 算法最大迭代次数
gain = 0.2; % CLEAN 环路增益 (0 < gain < 1)
threshold = 0.05; % 停止阈值（峰值低于最大值的比例）

%% 2. 生成仿真信号（目标 + 杂波 + 噪声）
n = 0:N-1;
% 目标信号（我们想保留的弱信号）
target_signal = exp(1i*2*pi*f0*n);
% 强杂波（需要抑制的干扰）
clutter_signal = zeros(1, N);
for k = 1:length(f_clutter)
    clutter_signal = clutter_signal + A_clutter(k)*exp(1i*2*pi*f_clutter(k)*n);
end
% 加性高斯白噪声
noise_power = 10^(-SNR_dB/10);
noise = sqrt(noise_power/2)*(randn(1,N)+1i*randn(1,N));
% 合成接收信号
received_signal = target_signal + clutter_signal + noise;

%% 3. CLEAN 算法核心实现
function[clean_signal, components, residue] = clean_algorithm(signal, K, gain, threshold)
    % signal: 输入信号（时域或频域，这里假设为时域）
    % K: 最大迭代次数
    % gain: 环路增益（控制每次减去多少比例的分量）
    % threshold: 停止阈值（当最大峰值小于初始最大峰值的 threshold 倍时停止）
    N = length(signal);
    residue = signal(:); % 初始残差 = 原始信号
    components = []; % 存储找到的分量 [位置，幅度，相位]
    
    % 计算初始信号频谱
    spec_original = fft(signal);
    max_peak_original = max(abs(spec_original));
    threshold_value = threshold * max_peak_original;
    
    for iter = 1:K
        % 3.1 计算当前残差信号的频谱
        spec_residue = fft(residue);
        % 3.2 找到频谱中的最大峰值及其位置
        [peak_value, peak_idx] = max(abs(spec_residue));
        % 3.3 检查是否达到停止条件
        if peak_value < threshold_value
            fprintf('迭代 %d: 峰值 (%.4f) 低于阈值 (%.4f)，停止。\n', ... iter, peak_value, threshold_value);
            break;
        end
        % 3.4 提取该频率分量的完整信息（幅度和相位）
        peak_phase = angle(spec_residue(peak_idx));
        peak_amplitude = abs(spec_residue(peak_idx))/ N; % 从 FFT 幅度转换为实际幅度
        % 3.5 存储找到的分量信息
        % 位置转换为归一化频率 (0~1)
        freq = mod(peak_idx-1, N);
        if freq > N/2
            freq = freq - N; % 转换为负频率表示
        end
        freq_norm = freq / N;
        components = [components; freq_norm, peak_amplitude, peak_phase, iter];
        % 3.6 从残差中减去该分量（乘以增益，防止过减）
        % 重建该分量的时域信号
        t = (0:N-1)';
        component_signal = peak_amplitude * exp(1i*(2*pi*freq_norm*t + peak_phase));
        % 更新残差
        residue = residue - gain * component_signal;
        % 显示迭代信息
        fprintf('迭代 %2d: 频率=%.4f, 幅度=%.4f, 残差峰值=%.4f\n', ... iter, freq_norm, peak_amplitude, peak_value);
    end
    % 4. 用找到的分量重建'干净'信号
    clean_signal = zeros(N,1);
    if ~isempty(components)
        t = (0:N-1)';
        for i=1:size(components,1)
            freq = components(i,1);
            amp = components(i,2);
            phase = components(i,3);
            clean_signal = clean_signal + amp * exp(1i*(2*pi*freq*t + phase));
        end
    end
end

%% 5. 运行 CLEAN 算法
tic;
[clean_signal, components, residue_signal] = clean_algorithm(received_signal.', K, gain, threshold);
time_elapsed = toc;
fprintf('CLEAN 算法完成，耗时 %.3f 秒，共找到 %d 个分量。\n\n', time_elapsed,size(components,1));

%% 6. 结果可视化
% 6.1 设置统一频率轴
freq_axis = (-N/2:N/2-1)/N;
spec_original = fftshift(abs(fft(received_signal)));
spec_clean = fftshift(abs(fft(clean_signal)));
spec_residue = fftshift(abs(fft(residue_signal)));
spec_target = fftshift(abs(fft(target_signal))); % 真实目标信号频谱

% 6.2 绘制频谱对比图
figure('Position',[100,100,1200,800]);
% 原始接收信号频谱
subplot(3,2,1);
plot(freq_axis, spec_original,'b','LineWidth',1.5);
xlabel('归一化频率');ylabel('幅度');title('(a) 原始接收信号频谱（含杂波与噪声）');xlim([-0.5,0.5]); grid on; hold on;
% 标记目标频率
plot([f0, f0],[0,max(spec_original)],'r--','LineWidth',1);
legend('接收信号','目标频率','Location','northwest');

% CLEAN 后信号频谱
subplot(3,2,2);
plot(freq_axis, spec_clean,'g','LineWidth',1.5);
xlabel('归一化频率');ylabel('幅度');title('(b) CLEAN 后信号频谱');xlim([-0.5,0.5]); grid on; hold on;
plot(freq_axis, spec_target,'r:','LineWidth',1.5);
legend('CLEAN 输出','真实目标','Location','northwest');

% 残差信号频谱
subplot(3,2,3);
plot(freq_axis, spec_residue,'m','LineWidth',1.5);
xlabel('归一化频率');ylabel('幅度');title('(c) 最终残差信号频谱');xlim([-0.5,0.5]); grid on;

% 频谱对比（重叠显示）
subplot(3,2,4);
plot(freq_axis, spec_original,'b:','LineWidth',1,'DisplayName','原始信号'); hold on;
plot(freq_axis, spec_clean,'g-','LineWidth',1.5,'DisplayName','CLEAN 后');
plot(freq_axis, spec_target,'r--','LineWidth',1.5,'DisplayName','真实目标');
xlabel('归一化频率');ylabel('幅度');title('(d) 频谱对比');xlim([-0.5,0.5]);legend('show'); grid on;

% 6.3 找到的分量列表
subplot(3,2,5);
if ~isempty(components)
    bar(components(:,1),components(:,2));
xlabel('归一化频率');ylabel('幅度');title('(e) CLEAN 找到的频率分量');xlim([-0.5,0.5]); grid on;
else
    text(0.5,0.5,'未找到分量','HorizontalAlignment','center'); axis off;
end

% 6.4 迭代收敛过程
subplot(3,2,6);
if ~isempty(components)
    plot(components(:,4),components(:,2),'ro-','LineWidth',1.5,'MarkerSize',8);
xlabel('迭代次数');ylabel('分量幅度');title('(f) 分量幅度随迭代变化'); grid on;
end
sgtitle(sprintf('CLEAN 算法杂波抑制仿真 (K=%d, gain=%.2f, threshold=%.2f)', K, gain, threshold),'FontSize',14);

%% 7. 性能评估
% 计算抑制前后目标频率附近的能量比
target_bin = round(f0*N)+1;
target_band = max(1, target_bin-3):min(N, target_bin+3);
% 原始信号中目标频带的能量
E_original_target = sum(spec_original(target_band).^2);
E_original_total = sum(spec_original.^2);
INR_original = 10*log10(E_original_target/(E_original_total - E_original_target + eps));
% CLEAN 后信号中目标频带的能量
E_clean_target = sum(spec_clean(target_band).^2);
E_clean_total = sum(spec_clean.^2);
INR_clean = 10*log10(E_clean_target/(E_clean_total - E_clean_target + eps));
fprintf('\n======= 性能评估 =======\n');
fprintf('原始信号目标干涉比 (INR): %.2f dB\n', INR_original);
fprintf('CLEAN 后目标干涉比 (INR): %.2f dB\n', INR_clean);
fprintf('INR 改善量：%.2f dB\n', INR_clean - INR_original);
% 计算与真实目标的相关系数
corr_coef = abs(corrcoef(target_signal, clean_signal.'));
fprintf('与真实目标信号的相关系数：%.4f\n',corr_coef(1,2));