当前位置：首页 > news >正文

利用小波变换对跳频信号进行参数估计

news 2025/9/18 9:41:47

跳频信号是一种扩频通信技术，其载波频率在多个频点之间按照伪随机序列跳变。小波变换因其良好的时频局部化特性，非常适合分析这类非平稳信号。

跳频信号模型

跳频信号可以表示为：

\[s(t) = A \exp\left[j(2\pi f_n t + \phi_n)\right], \quad t \in [t_n, t_{n+1}] \]

其中：

\(A\) 是信号幅度
\(f_n\) 是第\(n\)个 hop 的载波频率
\(\phi_n\) 是第\(n\)个 hop 的初始相位
\(t_n\) 是第\(n\)个 hop 的开始时间

小波变换基础

小波变换提供了信号在时间和频率域的联合表示：

\[W_s(a, b) = \frac{1}{\sqrt{|a|}} \int_{-\infty}^{\infty} s(t) \psi^*\left(\frac{t-b}{a}\right) dt \]

其中：

\(a\) 是尺度参数（与频率相关）
\(b\) 是平移参数（与时间相关）
\(\psi(t)\) 是小波母函数

MATLAB实现代码

%% 跳频信号参数估计 - 小波变换方法
clear; close all; clc;%% 1. 生成跳频信号
fs = 100e6;          % 采样率 100MHz
T_total = 100e-6;    % 总时长 100μs
t = 0:1/fs:T_total-1/fs;
N = length(t);% 跳频参数
hop_duration = 10e-6;           % 每个频率持续10μs
num_hops = floor(T_total/hop_duration);
hop_points = round(linspace(1, N, num_hops+1));% 频率集 (10-40MHz范围内的5个频率)
freq_set = [10, 15, 25, 30, 40] * 1e6;
freq_sequence = freq_set(randi(length(freq_set), 1, num_hops));% 生成跳频信号
signal = zeros(1, N);
phase = 0;
for i = 1:num_hopsstart_idx = hop_points(i);end_idx = hop_points(i+1)-1;segment_length = end_idx - start_idx + 1;% 当前段的信号signal(start_idx:end_idx) = exp(1j*(2*pi*freq_sequence(i)*t(start_idx:end_idx) + phase));% 更新相位（保持相位连续）phase = mod(2*pi*freq_sequence(i)*t(end_idx) + phase, 2*pi);
end% 添加噪声
SNR = 20; % 信噪比(dB)
signal_power = mean(abs(signal).^2);
noise_power = signal_power / (10^(SNR/10));
noise = sqrt(noise_power/2) * (randn(1, N) + 1j*randn(1, N));
signal_noisy = signal + noise;%% 2. 小波变换分析
% 选择小波
wavelet_name = 'cmor3-3'; % Complex Morlet小波，适合频率分析% 设置尺度范围（对应频率范围）
fmin = 5e6;  % 最小频率 5MHz
fmax = 50e6; % 最大频率 50MHz
scales = 1./(fmax:-0.5e6:fmin) * (centfrq(wavelet_name) * fs);% 执行连续小波变换
[cwt_coeffs, frequencies] = cwt(real(signal_noisy), scales, wavelet_name, 1/fs);
cwt_power = abs(cwt_coeffs).^2;%% 3. 参数估计
% 时频脊线提取（寻找每个时间点能量最大的频率）
[~, max_idx] = max(cwt_power);
estimated_freqs = frequencies(max_idx);% 跳变时刻检测 - 基于频率变化率
freq_diff = diff(estimated_freqs);
hop_threshold = 0.5 * max(abs(freq_diff)); % 自适应阈值
hop_indices = find(abs(freq_diff) > hop_threshold) + 1;% 估计跳频参数
estimated_hop_times = t(hop_indices);
estimated_durations = diff([0, estimated_hop_times, T_total]);
estimated_freq_seq = zeros(1, length(estimated_durations));for i = 1:length(estimated_durations)-1seg_start = hop_indices(i);seg_end = hop_indices(i+1)-1 if i < length(estimated_durations)-1 else N;seg_freqs = estimated_freqs(seg_start:seg_end);estimated_freq_seq(i) = mean(seg_freqs); % 使用均值作为该段的频率估计
end%% 4. 结果显示与评估
figure('Position', [100, 100, 1200, 800]);% 原始信号时域波形
subplot(3, 2, 1);
plot(t*1e6, real(signal));
title('原始跳频信号 (时域)');
xlabel('时间 (μs)');
ylabel('幅度');
grid on;% 含噪信号时域波形
subplot(3, 2, 2);
plot(t*1e6, real(signal_noisy));
title('含噪跳频信号 (时域)');
xlabel('时间 (μs)');
ylabel('幅度');
grid on;% 小波时频谱
subplot(3, 2, [3, 4]);
imagesc(t*1e6, frequencies/1e6, 10*log10(cwt_power));
axis xy;
xlabel('时间 (μs)');
ylabel('频率 (MHz)');
title('小波时频谱');
colorbar;
clim = get(gca, 'CLim');
set(gca, 'CLim', clim(2) + [-40 0]); % 调整颜色范围以突出主要成分% 提取的时频脊线
hold on;
plot(t*1e6, estimated_freqs/1e6, 'w-', 'LineWidth', 2);
plot(estimated_hop_times*1e6, estimated_freqs(hop_indices)/1e6, 'ro', 'MarkerSize', 8, 'LineWidth', 2);
hold off;% 频率序列比较
subplot(3, 2, 5);
stem(1:num_hops, freq_sequence/1e6, 'b', 'filled', 'LineWidth', 2);
hold on;
stem(1:length(estimated_freq_seq), estimated_freq_seq/1e6, 'r--', 'LineWidth', 1.5);
title('频率序列估计');
xlabel('Hop序号');
ylabel('频率 (MHz)');
legend('真实频率', '估计频率');
grid on;% 持续时间估计
subplot(3, 2, 6);
true_durations = hop_duration * ones(1, num_hops);
stem(1:num_hops, true_durations*1e6, 'b', 'filled', 'LineWidth', 2);
hold on;
stem(1:length(estimated_durations), estimated_durations(1:end-1)*1e6, 'r--', 'LineWidth', 1.5);
title('Hop持续时间估计');
xlabel('Hop序号');
ylabel('持续时间 (μs)');
legend('真实持续时间', '估计持续时间');
grid on;%% 5. 性能评估
fprintf('=== 跳频信号参数估计性能评估 ===\n');% 频率估计误差
freq_error = zeros(1, min(num_hops, length(estimated_freq_seq)));
for i = 1:length(freq_error)freq_error(i) = abs(estimated_freq_seq(i) - freq_sequence(i));
end
fprintf('频率估计平均误差: %.2f kHz\n', mean(freq_error)/1e3);
fprintf('频率估计最大误差: %.2f kHz\n', max(freq_error)/1e3);% 跳变时刻检测性能
if length(hop_indices) == num_hops-1fprintf('成功检测到所有 %d 个跳变时刻\n', num_hops-1);
elsefprintf('检测到 %d 个跳变时刻，实际应有 %d 个\n', length(hop_indices), num_hops-1);
end% 持续时间估计误差
duration_error = zeros(1, min(num_hops, length(estimated_durations)-1));
for i = 1:length(duration_error)duration_error(i) = abs(estimated_durations(i) - hop_duration);
end
fprintf('持续时间估计平均误差: %.2f μs\n', mean(duration_error)*1e6);
fprintf('持续时间估计最大误差: %.2f μs\n', max(duration_error)*1e6);%% 6. 改进的跳变时刻检测算法
% 基于小波系数模极大值的跳变检测
figure('Name', '改进的跳变检测算法');
[cwt_coeffs_real, ~] = cwt(real(signal_noisy), scales, wavelet_name, 1/fs);% 计算小波系数模的梯度
cwt_modulus = abs(cwt_coeffs_real);
modulus_gradient = abs(diff(cwt_modulus, 1, 2));% 寻找梯度极大值点作为跳变候选
[~, max_scale_idx] = max(cwt_modulus);
selected_gradient = zeros(1, size(modulus_gradient, 2));
for i = 1:size(modulus_gradient, 2)selected_gradient(i) = modulus_gradient(max_scale_idx(i), i);
end% 自适应阈值
gradient_threshold = 0.3 * max(selected_gradient);
improved_hop_indices = find(selected_gradient > gradient_threshold);% 去除过于接近的检测点
min_hop_interval = round(0.8 * hop_duration * fs);
improved_hop_indices = improved_hop_indices([true, diff(improved_hop_indices) > min_hop_interval]);subplot(2, 1, 1);
plot(t(1:end-1)*1e6, selected_gradient);
hold on;
plot(t(improved_hop_indices)*1e6, selected_gradient(improved_hop_indices), 'ro');
line([t(1), t(end)]*1e6, [gradient_threshold, gradient_threshold], 'Color', 'r', 'LineStyle', '--');
title('小波系数模梯度');
xlabel('时间 (μs)');
ylabel('梯度幅度');
legend('梯度', '检测到的跳变', '阈值');
grid on;subplot(2, 1, 2);
imagesc(t*1e6, frequencies/1e6, 10*log10(cwt_power));
axis xy;
hold on;
plot(t(improved_hop_indices)*1e6, estimated_freqs(improved_hop_indices)/1e6, 'wo', 'MarkerSize', 10, 'LineWidth', 2);
title('改进的跳变检测结果');
xlabel('时间 (μs)');
ylabel('频率 (MHz)');
colorbar;

算法原理与关键技术

1. 小波基函数选择

对于跳频信号分析，推荐使用Complex Morlet小波（'cmor'），因为：

具有良好的时频局部化特性
提供复数输出，可以同时获取幅度和相位信息
在频率估计方面有优异性能

2. 尺度选择策略

尺度参数与小波中心频率和实际频率的关系：

\[f = \frac{f_c}{a \cdot \Delta t} \]

其中：

\(f\) 是实际频率
\(f_c\) 是小波中心频率
\(a\) 是尺度参数
\(\Delta t\) 是采样间隔

3. 时频脊线提取

时频脊线提取是小波变换参数估计的核心：

% 方法1：直接取最大能量点
[~, max_idx] = max(cwt_power);
estimated_freqs = frequencies(max_idx);% 方法2：基于贪心算法的脊线跟踪
ridge = zeros(1, size(cwt_power, 2));
ridge(end) = find(cwt_power(:, end) == max(cwt_power(:, end)), 1);
for i = size(cwt_power, 2)-1:-1:1neighbors = max(1, ridge(i+1)-1):min(size(cwt_power, 1), ridge(i+1)+1);[~, idx] = max(cwt_power(neighbors, i));ridge(i) = neighbors(idx);
end
estimated_freqs = frequencies(ridge);

4. 跳变时刻检测

跳变时刻检测的几种方法：

频率差分法：检测频率变化的剧烈点
小波系数模极大值法：利用小波变换对奇异点的检测能力
多尺度综合法：在不同尺度上检测跳变点，然后综合判断

5. 参数估计算法优化

% 频率估计优化 - 基于相位信息
function refined_freq = refine_frequency_estimation(signal_segment, initial_est, fs)% 使用相位差分法提高频率估计精度phase = angle(signal_segment);phase_diff = diff(unwrap(phase));freq_estimates = phase_diff * fs / (2*pi);% 去除异常值freq_estimates = freq_estimates(abs(freq_estimates - initial_est) < initial_est/2);refined_freq = mean(freq_estimates);
end% 跳变时刻精确定位
function exact_hop_time = refine_hop_time(signal, rough_hop_idx, window_size, fs)% 在粗略估计的跳变点附近进行精细搜索search_start = max(1, rough_hop_idx - window_size);search_end = min(length(signal), rough_hop_idx + window_size);search_segment = signal(search_start:search_end);% 使用小波变换的模极大值精确定位[cwt_coeffs, ~] = cwt(real(search_segment), 1, 'gaus2', 1/fs);[~, exact_idx] = max(abs(cwt_coeffs));exact_hop_time = search_start + exact_idx - 1;
end

参考代码利用小波变换对跳频信号进行参数估计 www.youwenfan.com/contentcnh/55252.html