Fourier decomposition method

  • l5_933334
    了解作者
  • 8.3KB
    文件大小
  • zip
    文件格式
  • 0
    收藏次数
  • VIP专享
    资源类型
  • 0
    下载次数
  • 2022-05-23 10:56
    上传日期
FDM方法的MATLAB程序,适用于分解非线性非平稳信号。自己试过,可以使用的。
MatlabCodeOfTheFourierDecompositionMethodFDM.zip
  • MatlabCodeOfFDM
  • FMD_examples.m
    6.6KB
  • nspplote.m
    7.4KB
  • FMD_Low2High_High2LowSacnning.m
    14KB
内容介绍
%function [FIBFsLowToHigh FIBFsHighToLow EnergLekage]=FMD_Low2High_High2LowSacnning(xt,Fs,t) % % Calling sequence- % [FIBFsLowToHigh FIBFsHighToLow EnergLekage]=FMD_Low2High_High2LowSacnning(xt,Fs,t)% % All arguments are required % % Input- % xt - 1-D time series % Fs - Sampling Frequency % t - true time % % Output- % FIBFsLowToHigh - FIBFs from Low to high frequency scan % FIBFsHighToLow - FIBFs from High to low frequency scan % EnergLekage - energy leakage % Two, Time-frequency-Energy plot % % References: % Singh P, Joshi S.D., Patney R.K., Saha K., % The Fourier Decomposition Method for nonlinear and nonstationary time % series analysis, arXiv:1503.06675 [stat.ME]. % % code writer: Pushpendra Singh, June, 2014. % code writer: % % This is highly non optimzed initial vesrion of FDM, Matlab code. It can be optimzed. function [xt_recov_IMFsLowToHigh xt_recov_IMFsHighToLow EnergyLekage] = FMD_Low2High_High2LowSacnning(xt,Fs,t) global PS_PhaseUnwrap; global CentralDiff; CentralDiff=1; % central finite difference, e.g. for delta function PS_PhaseUnwrap=0; threshold=-0*(10)^(-1); % idealy should be zero. L=length(xt); if 1 % adding zero at last, can create discontinuity in signal. e.g. org emd example if rem(L,2) == 1 % odd L=L+1; % make it even for faster FFT %xt=[xt 0]; % can create discontinuity xt=[xt xt(end)]; % repeat last data end end NFFT=L; N=(NFFT); Xk=fft(xt,NFFT)/L; % k=1, Xk, is real, (2 to N/2) are complex, (N/2+1) is real, (N/2+2 to N) are complex conjugate of (N/2 to 2) nTotalHormonics=N/2; %if 1 % this is much better, seen by examples 1. %% xt_recov_IMFsLowToHigh=zeros(1,N)'; IFOfSignalIF=zeros(1,N)'; xt_AnalyticFIBF=zeros(1,N)'; xt_recov_IMFsLowToHigh(:,1)=(ones(1,N))*Xk(1); % first IMF DC init=2; p=2; while(init<=N/2) [kk xt_recov_FIBF xt_recov_AnalyticFIBF IFOfSignal]=getIMFsScanAllLowToHigh(Xk,Fs,init,nTotalHormonics,threshold); %% init=kk+1 xt_recov_IMFsLowToHigh(:,p)=xt_recov_FIBF; IFOfSignalIF(:,p-1)=IFOfSignal; xt_AnalyticFIBF(:,p-1)=xt_recov_AnalyticFIBF; p=p+1; end %subplot(2,1,1) sp_PlotTF(xt_AnalyticFIBF,IFOfSignalIF,t,Fs,Fs/2); title('FDM: Time-Frequncy-Energy estimate of FIBFs (LTH-FS)','FontSize',16,'FontName','Times'); nn=(1:N); tttmp=Xk((N/2)+1).*cos(pi*(nn-1)); % N/2+1 part of FFT, last part of DFT xt_recov_IMFsLowToHigh(:,end+1)=tttmp'; % N/2+1 part of FFT, last part of DFT %end %if 1 % this is much better, seen by examples 1 and 3. xt_recov_IMFsHighToLow=zeros(1,N)'; IFOfSignalIF=zeros(1,N)'; xt_AnalyticFIBF=zeros(1,N)'; mm=1; %tttmp=Xk((N/2)+1).*cos(pi*(1:N-1)); xt_recov_IMFsHighToLow(:,mm)=tttmp'; % N/2+1 part of FFT, last part of DFT final=N/2; while(final>=2) mm=mm+1; [kk xt_recov_FIBF xt_recov_AnalyticFIBF IFOfSignal]=getIMFsScanAllHighToLow(Xk,Fs,final,threshold); %% final=kk-1 xt_recov_IMFsHighToLow(:,mm)=xt_recov_FIBF; xt_AnalyticFIBF(:,mm-1)=xt_recov_AnalyticFIBF; IFOfSignalIF(:,mm-1)=IFOfSignal; end %subplot(2,1,2) %figure sp_PlotTF(xt_AnalyticFIBF,IFOfSignalIF,t,Fs,Fs/2); title('FDM: Time-Frequncy-Energy estimate of FIBFs (HTL-FS)','FontSize',16,'FontName','Times'); %% residue xt_recov_IMFsHighToLow(:,mm+1)=(ones(1,N))*Xk(1); % first IMF DC %end %Calculate Energy Leakge FMF_energy=0; for np=1:length(xt_recov_IMFsHighToLow(1,:)) FMF_energy=FMF_energy+sum(xt_recov_IMFsHighToLow(:,np).*xt_recov_IMFsHighToLow(:,np)); end SignalEnergy=sum(xt.*xt); EnergyLekage=SignalEnergy-FMF_energy; tt=1; %%% function [kk xt_recov_FIBF xt_recov_AnalyticFIBF IFOfSignal]=getIMFsScanAllLowToHigh(Xk,Fs,init,final,threshold) global PS_PhaseUnwrap; global CentralDiff; xt_recov_FIBF=0; xt_recov_AnalyticFIBF=0; N=length(Xk); %EndPoints=fix(length(Xk)*1/100); EndPoints=3; %threshold=-(10)^(-3); NbOfNegIF=(length(Xk)*1/100); % 1% is ok Analytic_zt=0; n=1:length(Xk); col=1; PositiveIMF=1; for kk=init:final Analytic_zt=Analytic_zt+(Xk(kk)).*exp((1i*2*pi*(kk-1)*(n-1)/N)); if PS_PhaseUnwrap IMFs_phase=sp_unwrap(angle(Analytic_zt)); else IMFs_phase=unwrap(angle(Analytic_zt)); end tmp1=(IMFs_phase); if CentralDiff %if 1 % central diff is better tmp=((tmp1(3:end)-tmp1(1:end-2))/2)*(Fs/(2*pi)); tmp=[tmp(1) tmp tmp(end)]; %else % five point diff %tmp=((8*tmp1(4:end-1)-8*tmp1(2:end-3)+tmp1(1:end-4)-tmp1(5:end))/12)*(Fs/(2*pi)); %tmp=[tmp(1) tmp(1) tmp tmp(end) tmp(end)]; %three point %tmp=((-3*tmp1(1:end-2)+4*tmp1(2:end-1)-tmp1(3:end))/2)*(Fs/(2*pi)); %tmp=[tmp tmp(end) tmp(end)]; %end else tmp=diff(tmp1)*(Fs/(2*pi)); % better derivative, e.g. use a wtooth & square waves as example tmp=[tmp tmp(end)]; %tmp=[tmp(2) tmp(2:end-1) tmp(end-1) tmp(end-1)]; % First & last value may be negative, so 2:end-1 end %if((sum(tmp(EndPoints:end-EndPoints)<0)>NbOfNegIF)) if((min(tmp(EndPoints:end-(EndPoints-1)))<threshold) && (max(xt_recov_FIBF)~=0) && (PositiveIMF==1) && (kk<=final)) % check -ve IF PositiveIMF=0; %break; recordIMF(:,col)=xt_recov_FIBF; recordAnalyticFIBFs(:,col)=xt_recov_AnalyticFIBF; recordIMF_IF(:,col)=xt_IF; % IF recordIMFkk(col)=kk-1; col=col+1; end if(min(tmp(EndPoints:end-(EndPoints-1)))>=threshold) % check +ve IF, First value may be negative so 2:end PositiveIMF=1; end xt_recov_FIBF=2*real(Analytic_zt); % *2 as taking only half of component xt_recov_AnalyticFIBF=2*Analytic_zt; xt_IF=tmp; % instantaneous Freq % First Value if(kk==init) recordIMF(:,col)=xt_recov_FIBF; recordAnalyticFIBFs(:,col)=xt_recov_AnalyticFIBF; recordIMF_IF(:,col)=xt_IF; % IF recordIMFkk(col)=kk; col=col+1; end %% for last value if((kk==final) && (PositiveIMF==1)) recordIMF(:,col)=xt_recov_FIBF; recordAnalyticFIBFs(:,col)=xt_recov_AnalyticFIBF; recordIMF_IF(:,col)=xt_IF; % IF recordIMFkk(col)=kk; end end getHighestValue=length(recordIMFkk); xt_recov_FIBF=recordIMF(:,getHighestValue); % return last value IMF xt_recov_AnalyticFIBF=recordAnalyticFIBFs(:,getHighestValue); IFOfSignal=recordIMF_IF(:,getHighestValue); kk=recordIMFkk(end); % return last value tt=1; function [kk xt_recov_FIBF xt_recov_AnalyticFIBF IFOfSignal]=getIMFsScanAllHighToLow(Xk,Fs,final,threshold) global PS_PhaseUnwrap; global CentralDiff; xt_recov_FIBF=0; xt_recov_AnalyticFIBF=0; N=
评论
    相关推荐
    • Matlab合集
      冈萨雷斯数字图像处理MATLAB版.中文版+数字图像处理第二版中文版(冈萨雷斯)+MATLAB-R2014a完全自学一本通+MATLAB R2016a完全自学一本通 素材文件+[模式识别与智能计算:MATLAB技术实现(第2版)].杨淑莹.扫描版
    • MATLAB教程
      MATLAB教程MATLAB教程MATLAB教程MATLAB教程MATLAB教程MATLAB教程
    • MATLAB
      MATLAB 该项目是在matlab上完成的,涉及创建和移动宇宙飞船和机器人。 太空飞船和机器人是使用Matlab中的简单几何形状创建的,并通过连续变换矩阵进行移动。 这个项目教我如何使用变换矩阵(旋转,平移等)的概念...
    • MATLAB基础
      一本学习matlab的一本好书
    • MATLAB编译器
      基于MATLAB 2018b版本介绍MATLAB编译器。介绍如何利用编译器将MATLAB代码编译为独立应用程序或组件,并在没有安装MATLAB的计算机上进行部署。
    • matlabruntime
      通过奇点容器运行您的matlab项目 可以在没有MATLAB的容器中运行matlab代码。 为此,我们首先需要通过Matlab编译器在本地构建相应的Matlab代码的独立应用程序,然后使用具有Matlab运行时( 的容器)运行该应用程序 )...
    • matlab实现
      matlab实现 matlab实现matlab实现matlab实现matlab实现
    • matlab 教程
      matlab 信号处理资料,里面包含信号处理pdf文档,一些杂乱的程序和命令等
    • matlab教程
      matlab教程,PPT格式,包含 matlab基本知识、matlab入门、matlab作图、线性规划、无约束优化、非线性规划、统计工具箱、差值、微分方程等多项知识点,并且每个知识点独立成为PPT,内还含有matlab信号处理详解等文档...
    • matlab简介
      1.MATLAB 开发环境 1.1 MATLAB 的视窗环境 进入MATLAB之后,会看到一个视窗MATLAB Command Window称为指令视窗,它是你键入指令的地方同时 MATLAB也将计算结果显示在此。 1.2 简易计算 我们先从MATLAB的...