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.