论文Chirp hunting分析

---

1．论文概述

这篇文章的目的是将一个信号分解为多个线性调频信号叠加的形式，即

其中

为高斯包络的线性调频信号。

2. 分析

在文章的第四章节需要估计信号的调频率，通过仿真发现论文所提出的基于LAF函数的调频率估计方法，效果非常不稳定，估计误差差别较大。
仿真采用的信号为单频点线性调频信号

s_t=exp(-1j*pi*K*t.^2);

仿真了文章第四章LAF的积分公式

当调频率设置为K=-10000时，积分值随着角度alpha的变化关系为

所得到的积分值最大的角度即为所估计信号调频率对应角度，其值为 -9.2180e+03。而当K=-5000时，结果为

估计结果失效。

3. 代码

主函数：

clc;clear;close all
%% 参数设置
B=1000;
T=0.1;
K=-B/T/2;
fs=2200;
t=(0:1/fs:T);
%% 信号采样
s_t=exp(1j*pi*K*t.^2);
[Spec,Freq,t] = Chirplet_Transform(s_t,length(s_t),round(length(s_t)/10),fs,K/2);
figure;colormap(jet);
imagesc(t,Freq,abs(Spec));
%% Wigner 分布
[tfr,f] = tfrWV(s_t.',length(t),fs);
figure;
colormap(jet);
imagesc(t,f,abs(tfr));
%% Chirp 分解
chirp_hunting(s_t,t,K,round(length(s_t)/10),fs)

Chirplet 变换:

function [Spec,Freq,t] = Chirplet_Transform(Sig,fLevel,WinLen,SampFreq,alpha)
% This function is used to calculate the Chirplet transform of the signal.
% --------------------INPUT------------------%
% Sig    一维原始信号 为横向量
% fLevel 频谱长度
% WinLen 高斯窗函数长度
% SampFreq 信号采样频率
% alpha  信号调频率（chirp rate）
% --------------------OUTPUT------------------%
% Spec   Chirplet变换结果
% Freq   垂直频率轴
% t      水平时间轴
%--------------------------------------------------------%
% written by Guowei Tu, 28/07/2019 in SME,SJTU (Contact me via GuoweiTu@sjtu.edu.cn)
%--------------------------------------------------------%
% the Original Chirplet Transform Algorithm is introduced in [1];
% while the codes here is motivated by [2], which provides a new insight to the Chirplet Transform.
%--------------------------------------------------------%
% [1].Steve Mann, Simon Haykin, The Chirplet Transform: Physical Considerations, 
% IEEE TRANSACTIONS ON SIGNAL PROCESSING,VOL. 43, NO. 11, NOVEMBER 1995
% [2].Peng Z.K , Meng G., Lang Z.Q.,Chu F.L, Zhang W.M., Yang Y., Polynomial Chirplet Transform with Application to Instantaneous Frequency Estimation,
% IEEE Transactions on Measurement and Instrumentation 60(2011) 3222-3229

if (nargin < 1)
    error('At least one parameter required!');
end
if (nargin < 4)       % 默认采样率 1000Hz
    SampFreq = 1000;
end
if (nargin < 3)       % 默认窗函数长度64
    WinLen = 64;
end
if (nargin < 2)       % 默认频谱长度512
    fLevel = 512;
end

%% data preparation
SigLen = length(Sig); % 信号长度
t = (0:1:SigLen-1)*(1/SampFreq); % 计算时间轴
%Sig = hilbert(Sig); % 将信号转换为解析信号（复数形式，只存在正频率，并且和原始信号正频率值相同）
%% Frequency axis and its half-axis points
fLevel = ceil(fLevel/2) * 2+1;    % 将频谱长度转化为最近的奇数
Lf = (fLevel - 1)/2;    % 频谱长度的一半
%% Generate Gauss window functions
WinLen = ceil(WinLen/2) * 2+1;    % 将窗函数长度转化为最接近的奇数
WinFun = exp(-6* linspace(-1,1,WinLen).^2);    % 高斯窗函数, 固定时间宽度 [-1, 1], 数据点长度为WinLen
Lw = (WinLen - 1)/2;    % 窗函数长度的一半
%% CT spectrum result array (initialized to 0)
Spec = zeros(fLevel,SigLen); %用于存储Chirplet变换结果
%% Sliding window signal,data segmentation preparation
for iLoop = 1:SigLen
    % Determine the upper and lower limits of the left and right signal index subscripts (note that to prevent the edge width from exceeding the time domain, the retrieval error!)
    iLeft = min([iLoop-1, Lw, Lf]);
    iRight = min([SigLen-iLoop, Lw, Lf]); %左右两边的滑动窗口均取这些数的最小值，以防止数组越界
    iIndex = -iLeft:iRight;%滑动窗口索引值
    
    iIndex1 = iIndex + iLoop;   % 原始信号的索引
    iIndex2 = iIndex + Lw + 1;  % 窗函数的索引
    Index = iIndex + Lf +1;     % 频率轴索引
    
    R_Operator = exp(-1i * 2*pi*alpha * t(iIndex1).^2 / 2); % 频率旋转算子
    S_Operator = exp(1i * 2*pi*alpha * t(iLoop) * t(iIndex1)); % 频率移位算子
    
    Sig_Section = Sig(iIndex1).*R_Operator.*S_Operator; % 信号滑动窗片段乘以旋转和移位算子，正常的Chirplet变换还应乘以相位exp(-1i * 2*pi*alpha * t(iLoop).^2 / 2),
                                                        %这里给省略了，因为对于每个给定的iLoop,该相位项为常数项。
    Spec(Index, iLoop) = Sig_Section.* conj(WinFun(iIndex2));  % 乘以窗函数
end
%% STFT
% Spec = fft(Spec); % 最后进行STFT
% Spec = Spec*2/fLevel; % restores the true magnitude of FT，这是由于对解析信号进行的操作，幅度减小了一半
% Spec = Spec(1:(end-1)/2,:); % 截止到采样频率的一半，这里猜测程序是为实信号而设计的
% fk = 0:1:fLevel-1;
% fk = fk(1:(end-1)/2);
% Freq = linspace(0,0.5*SampFreq,length(fk)); % 计算频率轴

Spec = fftshift(fft(Spec),1); % 最后进行STFT
Freq = linspace(-0.5*SampFreq,0.5*SampFreq,fLevel); % 计算频率轴
end

Wigner-Ville分布：

function [tfr,f] = tfrWV(x,N,SampFreq)
%	计算一维离散时间信号x的 Wigner-Ville distribution
 
%	x 原始信号，为列向量
%	N 频率单元个数
%	TFR 输出的时频分布矩阵
%	f 输出频率值
%   SampFreq 信号采样频率
%	refer to F. Auger, May-August 1994, July 1995.

xrow=length(x);%计算信号长度
t=(1:xrow);%信号长度的索引值
tfr= zeros(N,xrow); %初始化时频分布矩阵 

for icol=1:xrow
 ti= t(icol); 
 taumax=min([ti-1,xrow-ti,round(N/2)-1]);%取左余量、右余量以及频谱个数一半中的最小值
 tau=-taumax:taumax; 
 indices= rem(N+tau,N)+1;%使得时频矩阵行的索引永远为正
 tfr(indices,icol) = x(ti+tau) .* conj(x(ti-tau));
 tau=round(N/2); 
 if (ti<=xrow-tau)&&(ti>=tau+1)%当N=xrow时，该条件不成立
  tfr(tau+1,icol) = 0.5 * (x(ti+tau) * conj(x(ti-tau))  + ...
                           x(ti-tau) * conj(x(ti+tau))) ;
 end
end

tfr= fft(tfr); 
tfr=real(tfr); 
f=(0.5*(0:N-1)/N)'*SampFreq;
end

LAF函数：

function [alpha_3] = LAF(WV,I_f,I_t,r,a)
%LAF 用于实现论文Chirp hunting 中利用LAF(local ambiguity function)中
%求解调频率中的积分公式
%   WV 信号的魏格纳分布
%   I_f,I_t 待求时间频率点位置
%   a 角度变量
%   r 积分长度
%   alpha_3 积分公式结果
[x,y]=size(WV);
N_r=20; %取积分的个数
Q=zeros(length(a),N_r);%LAF函数
rr=linspace(-r,r,N_r);
delt_rr=rr(2)-rr(1);

for i=1:length(a)
 Index_t=max(min(rr*cos(a(i))+I_t),1):(delt_rr*cos(a(i))):min(max(rr*cos(a(i))+I_t),y);%设置索引用于防止数组越界
 Index_f= rem((Index_t-I_t)*tan(a(i))+I_f,x);
 n=length(Index_t);
 for j=1:n
    Q(i,j)=interp2(WV,Index_f(j),Index_t(j))*interp2(WV,Index_f(n-j+1),Index_t(n-j+1));
 end
end
alpha_3=sum(Q,2);

end

参考

J. C. O’Neill and P. Flandrin, “Chirp hunting,” Proceedings of the IEEE-SP International Symposium on Time-Frequency and Time-Scale Analysis (Cat. No.98TH8380), Pittsburgh, PA, USA, 1998, pp. 425-428, doi: 10.1109/TFSA.1998.721452.