# 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的...