?? tfar.m
字號:
function [AmlEst, B0lEst, EstInfo]= tfar(x, MMAX, LMAX, crit)% function [AmlEst, B0lEst, EstInfo]= tfar(x[, MMAX, LMAX, crit])% This file is part of the TFPM toolbox v1.0 (c)% michael.jachan@tuwien.ac.at and underlies the GPL.% % An automatic procedure to fit a TFAR(M, L) model to a non- % stationary signal x. The maximum search space defaults to % MMAX= min(8, floor(N/16)), % LMAX= min(4, floor(N/16)), % The order estimation IC <crit> defaults to the AIC. % EstInfo= {SigHat, InfoCr, NrPara, PENALTY, Ax};%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%clear;tfpm;N = 256;MAR = 3;LAR = 2;MMA = 0;LMA = LAR;re_im= 'i';mo_no= 'n';tfpm_file_gen;%-------------alpha= 1/2;beta = 1/2;x= tfarma_gen(randn(N, 1), Aml, Bml, beta);MMAX= 2*MAR;LMAX= 2*LAR;crit= 'AIC';P0= tfarma_wvsp(Aml, Bml, N, alpha);figure(1);tf_show(P0);drawnow%nargin= 1;% Some call examples [AmlEst, B0lEst, EstInfo]= tfar(x, MMAX, LMAX, '1.5');[AmlEst, B0lEst, EstInfo]= tfar(x, MMAX, LMAX);[AmlEst, B0lEst, EstInfo]= tfar(x, MMAX);[AmlEst, B0lEst]= tfar(x);%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Some constants, etc%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%alpha = 1/2; % TF shift parameterbeta = 1/2;N = length(x);lambda= .9800; % maximum root magnituderho = 1-log(12);% for the MDL%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% The nargin sory%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%if(nargin<=3) crit= 'AIC';end;if(nargin<=2) LMAX= min(8, floor(N/16));end;if(nargin<=1) MMAX= min(8, floor(N/16)); LMAX= min(4, floor(N/16));end;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% The penalty%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%switch(lower(crit)) case 'mdl' PENALTY= log(N+1)+rho; case 'aic' PENALTY= 2; otherwise PENALTY= str2num(crit);end;MMAX= min(MMAX, floor(N/16));LMAX= min(LMAX, floor(N/16));SigHat= exp(99)*ones(MMAX+1, LMAX+1);NrPara= zeros(MMAX+1, LMAX+1);InfoCr= NrPara;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Nonparametric signal statistics%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Ax= fft(corr_est(x, x, MMAX, 1/2));Ax= fft(corr_est(x, x, -1, 1/2));Ax= [Ax(N/2+1:N, :); Ax(1:N/2, :)];Rx= nm_to_nk(ml_to_nm(Ax));%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% TFAR Identification %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%sprintf('Searching up to TFAR(%d, %d) using the "%s"', MMAX, LMAX, crit)AR= {};for L= 0:LMAX warning off MATLAB:divideByZero Mmax= min(max(MMAX, floor(N/16/L)), MMAX); warning on MATLAB:divideByZero%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%if(0)% TAPER DESIGN (NOT ACTIVE)%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% [Psi, Mask, v2]= tf_multiwin(N, 3*Mmax, 4*L, 0, 2, 1); v2 subplot(1, N/4+1, L+1) tf_show(Psi);drawnow%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% [AAR, BAR]= tfar_est_tfywu(Ax(N/2+1-3*L:N/2+1+3*L, N/2+1-Mmax:N/2+1+Mmax), N); for M= 1:Mmax B= BAR(:, :, M);%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%if(0)% VARIANCE COMPUTATION BY INVERSE FILTERING%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% [Aml1, lambdamax, mmAR, nrFDIR]= param_stabilize(AAR(:, 1:M+1, M), N, lambda); [M L mmAR]; AR{M+1, L+1} = {Aml1 B}; eprime= tfarma_inv(x, Aml1, B/param_get(B, 0, 0));% figure(99);subplot(2, 1, 1);plot(real(eprime));subplot(2, 1, 2);plot(imag(eprime));title(sprintf('TFAR(%d, %d); %d', M, L, mmAR));drawnow SigHat(M+1, L+1)= eprime'*eprime/N; NrIter(M+1, L+1)= mmAR; EstInn(M+1, L+1, :)= eprime;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%if(0)% VARIANCE COMPUTATION BY INNER PRODUCT%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% AR{M+1, L+1} = {AAR(:, 1:M+1, M) B}; Lalpha= tfarma_weyl(AAR(:, 1:M+1, M), B/param_get(B, 0, 0), N, alpha); SigHat(M+1, L+1)= real(sum(sum(conj(Rx)./abs(Lalpha)^2)))/N;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% end;end;%for L= 0:LMAXNrPara(:, :)= (1:MMAX+1)'*(2*(0:LMAX)+1) - 1;InfoCr(:, :)= log(SigHat(:, :)) + PENALTY*NrPara(:, :)/N;[MAREst, LAREst] = find(InfoCr(:, :)==min(min(InfoCr(:, :))));MAREst= MAREst-1;LAREst= LAREst-1;AmlEst= param_stabilize(AR{MAREst+1, LAREst+1}{1}, N, lambda);B0lEst= AR{MAREst+1, LAREst+1}{2};%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Further Detailled Info%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%if(nargout==3) EstInfo= {SigHat, InfoCr, NrPara, PENALTY, Ax};end;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%clear;tfpm;N = 256;MAR = 3;LAR = 2;MMA = 0;LMA = LAR;re_im= 'i';mo_no= 'n';tfpm_file_gen;%-------------alpha= 1/2;beta = 1/2;MMAX= 6;LMAX= 6;MM= 20;STATUS= {};for mm= 1:MM mm x= tfarma_gen(randn(N, 1), Aml, Bml); [AmlEst, B0lEst, EstInfo]= tfar(x, MMAX, LMAX, 'MDL'); STATUS{1, mm}= EstInfo; [AmlEst, B0lEst, EstInfo]= tfar(x, MMAX, LMAX, 'AIC'); STATUS{2, mm}= EstInfo;end;filename= sprintf('tfar_%04d%d%d%s%s', N, MAR, LAR, re_im, mo_no)save(filename, 'STATUS')%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%clear;tfpm;N = 256;MAR = 3;LAR = 2;MMA = 1;LMA = LAR;re_im= 'i';mo_no= 'n';tfpm_file_gen;%-------------alpha= 1/2;beta = 1/2;mm= 1;MM= 20;filename= sprintf('tfarma20_%04d%d%d%d%d%s%s', N, MAR, LAR, MMA, LMA, re_im, mo_no)load(filename)TFAR= 0;TFMA= 0;TFAA= 0;MAR= 0;LAR= 0;MMA= 0;LMA= 0;MAA= 0;LAA= 0;for mm= 1:MM InfoCr= STATUS{2, 2, mm}{2}; [MAREst, LAREst] = find(InfoCr(:, :, 1)==min(min(InfoCr(:, :, 1)))); MAREst= MAREst-1;LAREst= LAREst-1; [MMAEst, LMAEst] = find(InfoCr(:, :, 2)==min(min(InfoCr(:, :, 2)))); MMAEst= MMAEst-1;LMAEst= LMAEst-1; [MAAEst, LAAEst] = find(InfoCr(:, :, 3)==min(min(InfoCr(:, :, 3)))); MAAEst= MAAEst-1;LAAEst= LAAEst-1; if(InfoCr(MAREst+1, LAREst+1, 1)<InfoCr(MMAEst+1, LMAEst+1, 2) & ... InfoCr(MAREst+1, LAREst+1, 1)<InfoCr(MAAEst+1, LAAEst+1, 3)) sprintf('TFAR(%d, %d)', MAREst, LAREst) MAR= MAR+MAREst;LAR= LAR+LAREst; TFAR= TFAR+1; else if(InfoCr(MMAEst+1, LMAEst+1, 2)<InfoCr(MAAEst+1, LAAEst+1, 3)) sprintf('TFMA(%d, %d)', MMAEst, LMAEst) MMA= MMA+MMAEst;LMA= LMA+LMAEst; TFMA= TFMA+1; else sprintf('TFARMA(%d, %d)', MAAEst, LAAEst) MAA= MAA+MAAEst;LAA= LAA+LAAEst; TFAA= TFAA+1; end; end;end;[100*TFAR/MM MAR/TFAR LAR/TFAR][100*TFMA/MM MMA/TFMA LMA/TFMA][100*TFAA/MM MAA/TFAA LAA/TFAA]%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%?? 快捷鍵說明
復制代碼
Ctrl + C
搜索代碼
Ctrl + F
全屏模式
F11
切換主題
Ctrl + Shift + D
顯示快捷鍵
?
增大字號
Ctrl + =
減小字號
Ctrl + -