CPF实现记录

1. 工作步骤

• 弄清楚WVD实现源码，以及横坐标，纵坐标的变换，阈值等等。
• 用matlab时频工具箱实现信号的WVD变换。（具体参数下文给出 ，以下都是）
• 找到matlab时频工具箱实现WVD变换的源码。
• 然后根据找到的源码，参考源码，写出信号的模糊函数分析AF，
• 参考WVD源码实现信号的CPF。

2. WVD实现需要注意事项

1. 信号得是列信号：

...
t=0:1/fs:T-1/fs;
t=t';
....
X=A1*exp(1i*(a_10+a_11*t+a_12*t.^2));
[s_lfm_wv,t,f]=tfrwv(X);
contour(t/fs,f*fs,abs(s_lfm_wv));

得到得图形：

如果使用列信号：

...
t=0:1/fs:T-1/fs;
....
X=A1*exp(1i*(a_10+a_11*t+a_12*t.^2));
[s_lfm_wv,t,f]=tfrwv(X');
contour(t/fs,f*fs,abs(s_lfm_wv));

得到得图形为：

从图形来看，和我们要得最终结果是相反得，斜率是负得，所以需要把纵坐标显示为250-0；

我这里就直接使用第一种情况t=t’了不翻转Y纵；

3. WVD 实现

open tfrwv可以看tfrwv.m源码：
根据源码得到多分量的LFM信WVD变换为：

核心代码：

%=================开始wvd变换代码==================
t_wvd=1:N;
[trow,tcol]=size(t_wvd);
[xrow,xcol]=size(X);
tfr=zeros(N,tcol);
for icol=1:tcol,
    t_wvdi=t_wvd(icol);
    taumax=min([t_wvdi-1,xrow-t_wvdi,round(N/2)-1]);
    tau=(-taumax:taumax);
    indices=rem(N+tau,N)+1;
    tfr(indices,icol)=X(t_wvdi+tau,1).*conj(X(t_wvdi-tau,xcol));
    tau=round(N/2);
    if(t_wvdi<=xrow-tau)&(t_wvdi>=tau+1),
        tfr(tau+1,icol)=0.5*(X(t_wvdi+tau,1)*conj(X(t_wvdi-tau,xcol))+...
           X(t_wvdi-tau,1)*conj(X(t_wvdi+tau,xcol)));
    end 
end
tfr1=fft(tfr,nfft,1);
if(xcol==1),tfr1=real(tfr1);end
f=(0.5*(0:N-1)/N)';
figure;
contour(t_wvd/fs,f,abs(tfr1));
title(['WVD变换',' f1=',num2str(f1),' ',' k1=',num2str(k1),' f2= ',num2str(f2),' k2= ',num2str(k2)]);

4. AF实现

同样的思路利用open ambifunb查看ambifunb.m源码，这里分为窄带AF和宽带AF

核心代码：使用了时频分析工具箱函数ambifunb.m

[NAF,TAU,XI]=ambifunb(X);
figure;
contour(TAU/fs,XI,abs(NAF));
xlabel('时间t(s)');
ylabel('归一化频率');
title(['窄带AF变换',' f1=',num2str(f1),' ',' k1=',num2str(k1),' f2= ',num2str(f2),' k2= ',num2str(k2)]);

二、 STFT相关仿真