【非线性线性调频模式分解】具有交叉信号模式的色散信号的群延迟（色散曲线）估计和模式分离（Matlab代码实现）

 👨‍🎓个人主页：研学社的博客  

💥💥💞💞欢迎来到本博客❤️❤️💥💥

🏆博主优势：🌞🌞🌞博客内容尽量做到思维缜密，逻辑清晰，为了方便读者。

⛳️座右铭：行百里者，半于九十。

📋📋📋本文目录如下：🎁🎁🎁

目录

💥1 概述

📚2 运行结果

2.1 色散信号

2.2 脉冲信号

🎉3 参考文献

🌈4 Matlab代码实现

---

💥1 概述

文献来源：

变分模态分解（VMD）是近年来引入的一种自适应数据分析方法，在各个领域引起了广泛的关注。然而，VMD是基于信号模型窄带特性的假设而制定的。为了分析宽带非线性线性调频信号（NCS），我们提出了一种称为变分非线性线性调频模式分解（VNCMD）的替代方法。VNCMD的开发源于宽带NCS可以通过使用解调技术转换为窄带信号的事实。因此，我们的分解问题被表述为最优解调问题，通过乘法器的交替方向方法可以有效地求解。我们的方法可以看作是一个时频滤波器组，它同时提取所有信号模式。提供了一些模拟和真实的数据示例，展示了VNCMD在分析包含接近甚至交叉模式的NCS方面的有效性。

📚2 运行结果

2.1 色散信号

 

 

2.2 脉冲信号

部分代码：

clc
clear
close all
SampFreq=10000;
load impulse.mat
T = 0.08; % Time duration
Nt = length(impulse);% The number of samples in time domain
Nf = floor(Nt/2)+1;% The number of samples in frequency domain
f = (0 : Nf-1)/T; % frequency variables
t = (0 : Nt-1)/SampFreq;  % time variables

Sign = awgn(impulse,0,'measured');

figure
set(gcf,'Position',[20 100 640 500]);        
set(gcf,'Color','w'); 
plot(t,Sign,'linewidth',2);
xlabel('Time (s)','FontSize',24,'FontName','Times New Roman');
ylabel('Amplitude (AU)','FontSize',24,'FontName','Times New Roman');
set(gca,'FontSize',24)
set(gca,'linewidth',2);
set(gca,'FontName','Times New Roman')

fftspec = fft(Sign); % Obtain frequency-domain data of the signal
Dsn = fftspec(1:Nf); 
%% SFFT
Nfrebin = 1024;
window = 64;
figure
[Spec,f] = SFFT(Dsn(:),SampFreq,Nfrebin,window);
imagesc(t,f,abs(Spec));  
axis([0 max(t) 0 SampFreq/2]);
set(gcf,'Position',[20 100 320 250]);
set(gcf,'Color','w');    
xlabel('Time (s)','FontSize',12,'FontName','Times New Roman');
ylabel('Frequency (Hz)','FontSize',12,'FontName','Times New Roman');
set(gca,'YDir','normal')
set(gca,'FontSize',12);

%% Parameter setting
beta = 1e-8;  % this parameter can be smaller which will be helpful for the convergence, but it may cannot properly track fast varying GDs
alpha = 3e-7;  % if this parameter is larger, it will help the algorithm to find correct modes even the initial GDs are too rough. But it will introduce more noise and also may increase the interference between the signal modes
tol = 1e-8;
%% Mode 1 extraction
Envelope = abs(hilbert(Sign)); % envelope signal of the impulse signal
[~,tindex1] = max(Envelope);
peakenv1 = t(tindex1); % GD initialization by finding peak time of the envelope signal
iniGD1 = peakenv1*ones(1,length(Dsn)); % initial GD vector
[eGDest1 temp1] = GDMD(Dsn,T,iniGD1,alpha,beta,tol); % extract the signal mode
Desest1 = temp1(1,:,end);
DFs1 = [Desest1,conj(fliplr(Desest1(2:ceil(Nt/2))))]; eM1 = real(ifft(DFs1)); %Obtain the time-domain signal by inverse FFT; DFs1 denotes the bilateral spectrums
%% Mode 2 extraction
Dsresidue1 = Dsn - Desest1; % obtain the residual signal by removing the extracted component from the raw signal (in frequency domain)
residue1 = Sign - eM1; % Residual signal in time domain
reEnvelope1 = abs(hilbert(residue1)); % envelope signal of the residual signal
[~,tindex2] = max(reEnvelope1);
peakenv2 = t(tindex2); % GD initialization by finding peak time of the envelope signal
iniGD2 = peakenv2*ones(1,length(Dsn)); % initial GD vector
[eGDest2 temp2] = GDMD(Dsresidue1,T,iniGD2,alpha,beta,tol); % extract the signal mode
Desest2 = temp2(1,:,end);
DFs2 = [Desest2,conj(fliplr(Desest2(2:ceil(Nt/2))))]; eM2 = real(ifft(DFs2)); %Obtain the time-domain signal by inverse FFT;
%% Mode 3 extraction
Dsresidue2 = Dsresidue1 - Desest2; % obtain the residual signal by removing the extracted component from the raw signal (in frequency domain)
residue2 = residue1 - eM2; % Residual signal in time domain
reEnvelope2 = abs(hilbert(residue2)); % envelope signal of the residual signal
[~,tindex3] = max(reEnvelope2);
peakenv3 = t(tindex3); % GD initialization by finding peak time of the envelope signal
iniGD3 = peakenv3*ones(1,length(Dsn)); % initial GD vector
[eGDest3 temp3] = GDMD(Dsresidue2,T,iniGD3,alpha,beta,tol); % extract the signal mode
Desest3 = temp3(1,:,end);
DFs3 = [Desest3,conj(fliplr(Desest3(2:ceil(Nt/2))))]; eM3 = real(ifft(DFs3)); %Obtain the time-domain signal by inverse FFT;
%%

eSig = eM1 + eM2 + eM3; % Reconstructed signal modes

figure
set(gcf,'Position',[20 100 640 500]);        
set(gcf,'Color','w'); 
plot(t,impulse,'linewidth',2);  
hold on
plot(t,eSig,'r--','linewidth',2); 
xlabel('Time (s)','FontSize',24,'FontName','Times New Roman');
ylabel('Amplitude (AU)','FontSize',24,'FontName','Times New Roman');
set(gca,'FontSize',24)
set(gca,'linewidth',2);
legend('True','Estimated')

%%%%% In practice, the above procedure can be executed iteratively until the energy of the residual signal is smaller than a threshold

🎉3 参考文献

部分理论来源于网络，如有侵权请联系删除。

🌈4 Matlab代码实现