basic_wvd_hough_single_signal.m

来自「这是关于LFM信号利用WHT变换估计其参数的仿真代码。程序基于时频分析包」· M 代码 · 共 89 行

% 进行基本的WVD和Hough变换尝试
% 利用Hough变换实现了对LFM信号初始频率和信号斜率的准确估计
% 程序功能验证正确，这里只指对单个信号处理，而没有涉及处理多个LFM信号
clc;
clear all;
close all;
%产生线性调频信号
N=512;
Start_Freq = -0.2;
End_Freq = 0.2;
% 注意，利用这个函数产生LFM信号时
% 信号对应的时间是1s，采样点数是N点，因此，相当于采样频率为1/N
% Start_Freq和End_Freq都是归一化频率，产生的LFM信号的斜率G=End_Freq-Start_Freq
sig = fmlin(N, Start_Freq, End_Freq);

% 信号时域波形
figure;
plot(1:N,real(sig),'LineWidth',2);
xlabel('时间 t');
ylabel('幅值 A');
title('产生信号的时域波形');

%%
% 计算WVD分布
% 这里对原始的WVD变换函数tfrwv进行了修改，得到tfrwv_Polo函数，修改代码部分具体参考函数内部注释
% 函数返回结果f是坐标结果，代表的频率范围是-0.25～0.25. 
% 对于FFT变换，频率轴对应的频率范围是-0.5～0.5，对于WVD变换，则变为它的一半-0.25～0.25
% 如果设置的信号频率超过此范围，则会发生周期旋转
[tfr_wvd,t,f]=tfrwv_Polo(sig, (1:N), N);
% 显示WVD变换
figure;
subplot(2,1,1);
contour(t, f, tfr_wvd);
xlabel('时间 t');
ylabel('频率 f');
title('wigner变换的等高图');
grid on;
subplot(2,1,2);
mesh(t, f, tfr_wvd);
shading('flat')
xlabel('时间 t');
ylabel('频率 f');
title('wigner变换的三维图');
grid on;


%Hough变换以及显示
% 对原始htl进行了修改，得到了htl_Polo函数，具体修改部分参考原始代码注释部分
[WH_wvd,rho,theta]=htl_Polo(tfr_wvd,N,N);
% 查找最大值点和对应的线
WH_wvd = ceil(WH_wvd);
P = houghpeaks(WH_wvd,1);
%显示Wigner-Hough分布等高线
figure;
subplot(2,1,1);
contour(theta, rho, WH_wvd);
hold on;
plot(theta(P(:,2)), rho(P(:,1)), 's', 'color', 'red');
xlabel('theta');
ylabel('rho');
title('wigner-hough变换的等高图');
grid on;
subplot(2,1,2);
mesh(theta, rho, WH_wvd);
shading('flat')
xlabel('theta');
ylabel('rho');
title('wigner-hough变换的三维图');
grid on;

%% 对线性调频信号进行估计
Signal_Accurate_G = End_Freq - Start_Freq;
% 这里的估计算法很重要，注意，这里考虑了一个1/2因子，原因是在于WVD变换时，利用的是FFT，
% 所得的频率坐标实际是原始频率值的两倍，Hough按照WVD图形得到那个图像对应的斜率结果
% 因此，根据那个WVD变换结果进行Hough变换检测得到的直线斜率和真实的LFM信号斜率之间的关系就是1/2的关系了。
Signal_Estimate_G = -1/(2*tan(theta(P(:,2))));
fprintf('实际信号斜率: %5.9f\n', Signal_Accurate_G);
fprintf('WHT检测得到的角度(弧度): %5.9f, 检测得到的信号斜率: %5.9f\n', theta(P(:,2)), Signal_Estimate_G);
fprintf('斜率估计误差: %5.9f\n\n', Signal_Accurate_G - Signal_Estimate_G);

Start_Freq_Accurate = Start_Freq;
% 这里的估算算法很重要，注意，这里考虑了一个平移。
% 因为，利用round(rho(P(:,1))/sin(theta(P(:,2))))只能得到t=N/2时刻的信号频率，而无法得到t=0时刻的信号频率
% 因此，这里估计起始频率的算法实际是利用t=N/2时刻的频率与斜率的估计对t=0时刻的频率进行估计。
% 具体算法就是斜率的1/2(因为N对应1，N/2对应就是1/2) 加上 利用round(rho(P(:,1))/sin(theta(P(:,2))))估计的t=N/2时刻的信号频率。
Start_Freq_Estimate = f(N/2-round(rho(P(:,1))/sin(theta(P(:,2))))) - 1/2*Signal_Estimate_G;
fprintf('信号起始频率: %5.9f\n', Start_Freq_Accurate);
fprintf('WHT检测得到的起始频率: %5.9f\n', Start_Freq_Estimate);
fprintf('初始频率估计误差: %5.9f\n\n', Start_Freq_Accurate - Start_Freq_Estimate );