?? tfarma_fit.m
字號:
function [MDL, BIC, AIC, GIC, MIC, mdl, bic, aic, gic, mic]= tfarma_fit(x, Mmax, Lmax, Psi, lambda)% function [MDL, BIC, AIC, GIC, mdl, bic, aic, gic]= tfarma_fit(x, Mmax, Lmax, Psi, lambda)% This file is part of the TFPM toolbox v1.0 (c)% michael.jachan@tuwien.ac.at and underlies the GPL.% % Searches TFARMA(M, L, M-1, L) models for M= 1:Mmax, L= 0:Lmax. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%clear;tfpm;N = 64;MAR = 2;LAR = 1;MMA = 1;LMA = 1;re_im= 'r';mo_no= 'n';tfpm_file_gen;%-------------alpha= 1/2;beta = 1/2;lambda= .9800;Mmax= max([2*MAR 2]);Lmax= max([2*LAR 2]);Psi= ones(N);%tf_multiwin(N, 3*Mmax, 4*Lmax, 0, 2, 1);x= tfarma_gen(randn(N, 1), Aml, Bml, beta);%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%alpha= 1/2;beta = 1/2;N= length(x);Ax= fft(corr_est(x, x, -1, alpha));Ax= [ Ax(N/2+1:N, :); Ax(1:N/2, :)].*conj(Psi);E= zeros(Mmax, Lmax+1);Eprime= zeros(Mmax, Lmax+1);for L= 0:Lmax for M= 1:Mmax AA= tfarma_est_tfywu(Ax(N/2-3*L+1:N/2+3*L+1, N/2+(M-1)-M+1:N/2+(M-1)+M+1), M-1, N);% Inverse TFAR Filtering [Aml1, lambdamax, mm]= param_stabilize(AA(:, 1:M+1, M), N, lambda); uM= tfarma_inv(x, Aml1, 1); Au= fft(corr_est(uM, uM, -1, alpha)); Au= [ Au(N/2+1:N, :); Au(1:N/2, :) ].*conj(Psi); [B, REG, PMIN]= tfma_est_cepsb(Au, M-1, L); [Bml1, lambdamax, mm]= param_stabilize(B, N, lambda); % Inverse TFMA Filtering eprime= tfarma_inv(uM, 1, Bml1/param_get(Bml1, 0, 0)); e = tfarma_inv(uM, 1, Bml1); E(M, L+1) = e'*e/N; Eprime(M, L+1)= eprime'*eprime/N; end;end;e= log(E);eprime= log(Eprime);rho= 1-log(12);MDL= eprime + (log(N+1)+rho) *(NA+NB+1/2)/N;BIC= eprime + log(N) *(NA+NB )/N;AIC= eprime + 2 *(NA+NB )/N;GIC= eprime + (2+rho) *(NA+NB )/N;MIC= eprime + (1+rho+log(N)/2)*(NA+NB )/N;mdl= e + (log(N+1)+rho) *(NA+NB+1/2)/N;bic= e + log(N) *(NA+NB )/N;aic= e + 2 *(NA+NB )/N;gic= e + (2+rho) *(NA+NB )/N;mic= e + (1+rho+log(N)/2)*(NA+NB )/N;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%clear;tfpm;for n= 8:8 N= 2^n; re_im= 'i'; mo_no= 'n'; alpha= 1/2; beta = 1/2; lambda= .9800; MM= 20; for M= 1:1 for L= 2:2 MAR = M; LAR = L; MMA= MAR-1; LMA = LAR; tfpm_file_gen;%------------------------- Mmax= max([2*M 2]); Lmax= max([2*L 2]); Psi= ones(N); [Psi, Mask, v2]= tf_multiwin(N, 2*Mmax+2, 2*Lmax+2, 5, 0, 0); v2 MDLhist= zeros(Mmax, Lmax+1); BIChist= zeros(Mmax, Lmax+1); AIChist= zeros(Mmax, Lmax+1); GIChist= zeros(Mmax, Lmax+1); MIChist= zeros(Mmax, Lmax+1); mdlhist= zeros(Mmax, Lmax+1); bichist= zeros(Mmax, Lmax+1); aichist= zeros(Mmax, Lmax+1); gichist= zeros(Mmax, Lmax+1); michist= zeros(Mmax, Lmax+1); for mm= 1:MM x= tfarma_gen(randn(N, 1), Aml, Bml/param_get(Bml, 0, 0), 1/2); [MDL, BIC, AIC, GIC, MIC, mdl, bic, aic, gic, mic]= tfarma_fit(x, Mmax, Lmax, Psi, lambda); [MMDL, LMDL] = find(MDL==min(min(MDL)));LMDL= LMDL-1; [MBIC, LBIC] = find(BIC==min(min(BIC)));LBIC= LBIC-1; [MAIC, LAIC] = find(AIC==min(min(AIC)));LAIC= LAIC-1; [MGIC, LGIC] = find(GIC==min(min(GIC)));LGIC= LGIC-1; [MMIC, LMIC] = find(MIC==min(min(MIC)));LMIC= LMIC-1; [Mmdl, Lmdl] = find(mdl==min(min(mdl)));Lmdl= Lmdl-1; [Mbic, Lbic] = find(bic==min(min(bic)));Lbic= Lbic-1; [Maic, Laic] = find(aic==min(min(aic)));Laic= Laic-1; [Mgic, Lgic] = find(gic==min(min(gic)));Lgic= Lgic-1; [Mmic, Lmic] = find(mic==min(min(mic)));Lmic= Lmic-1; [N MAR LAR mm MMDL, LMDL, MBIC, LBIC, MAIC, LAIC, MGIC, LGIC, MMIC, LMIC] [N MAR LAR mm Mmdl, Lmdl, Mbic, Lbic, Maic, Laic, Mgic, Lgic, Mmic, Lmic] MDLhist(MMDL, LMDL+1)= MDLhist(MMDL, LMDL+1) + 1; BIChist(MBIC, LBIC+1)= BIChist(MBIC, LBIC+1) + 1; AIChist(MAIC, LAIC+1)= AIChist(MAIC, LAIC+1) + 1; GIChist(MGIC, LGIC+1)= GIChist(MGIC, LGIC+1) + 1; MIChist(MMIC, LMIC+1)= MIChist(MMIC, LMIC+1) + 1; mdlhist(Mmdl, Lmdl+1)= mdlhist(Mmdl, Lmdl+1) + 1; bichist(Mbic, Lbic+1)= bichist(Mbic, Lbic+1) + 1; aichist(Maic, Laic+1)= aichist(Maic, Laic+1) + 1; gichist(Mgic, Lgic+1)= gichist(Mgic, Lgic+1) + 1; michist(Mmic, Lmic+1)= michist(Mmic, Lmic+1) + 1; figure(1);subplot(2, 2, 1);imagesc(MDLhist) set(gca, 'XTick', 1:Lmax+1) set(gca, 'XTickLabel', 0:Lmax) figure(1);subplot(2, 2, 2);imagesc(BIChist) set(gca, 'XTick', 1:Lmax+1) set(gca, 'XTickLabel', 0:Lmax) figure(1);subplot(2, 2, 3);imagesc(AIChist) set(gca, 'XTick', 1:Lmax+1) set(gca, 'XTickLabel', 0:Lmax) figure(1);subplot(2, 2, 4);imagesc(GIChist) set(gca, 'XTick', 1:Lmax+1) set(gca, 'XTickLabel', 0:Lmax) figure(2);subplot(2, 2, 1);imagesc(mdlhist) set(gca, 'XTick', 1:Lmax+1) set(gca, 'XTickLabel', 0:Lmax) figure(2);subplot(2, 2, 2);imagesc(bichist) set(gca, 'XTick', 1:Lmax+1) set(gca, 'XTickLabel', 0:Lmax) figure(2);subplot(2, 2, 3);imagesc(aichist) set(gca, 'XTick', 1:Lmax+1) set(gca, 'XTickLabel', 0:Lmax) figure(2);subplot(2, 2, 4);imagesc(gichist) set(gca, 'XTick', 1:Lmax+1) set(gca, 'XTickLabel', 0:Lmax) drawnow end; exp= 'ord'; resultname= sprintf('data/%04d/%s%d%d%s%s.mat', N, exp, M, L, re_im, mo_no) save(resultname, 'lambda', 'Mmax', 'Lmax', 'MM', 'MDLhist', 'BIChist', 'AIChist', 'GIChist', 'MIChist', 'mdlhist', 'bichist', 'aichist', 'gichist', 'michist') end; end;end;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%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;re_im= 'i';mo_no= 'n';alpha= 1/2;beta = 1/2;M= 3L= 2MSEn= [];meanMn= [];meanLn= [];for n= 6:9 N= 2^n; exp= 'ord'; resultname= sprintf('data/%04d/%s%d%d%s%s.mat', N, exp, M, L, re_im, mo_no) load(resultname) MDLM= sum(MDLhist.').';MDLL= sum(MDLhist); BICM= sum(BIChist.').';BICL= sum(BIChist); AICM= sum(AIChist.').';AICL= sum(AIChist); GICM= sum(GIChist.').';GICL= sum(GIChist); MICM= sum(MIChist.').';MICL= sum(MIChist); mdlM= sum(mdlhist.').';mdlL= sum(mdlhist); bicM= sum(bichist.').';bicL= sum(bichist); aicM= sum(aichist.').';aicL= sum(aichist); gicM= sum(gichist.').';gicL= sum(gichist); micM= sum(michist.').';micL= sum(michist);figure(n);clf;subplot(4, 2, 1);bar(MDLM, 1);axis([.5 Mmax+0.5 0 MM])subplot(4, 2, 2);bar(MDLL, 1);axis([.5 Lmax+1.5 0 MM])set(gca, 'XTick', 1:Lmax+1)set(gca, 'XTickLabel', 0:Lmax)subplot(4, 2, 3);bar(BICM, 1);axis([.5 Mmax+0.5 0 MM])subplot(4, 2, 4);bar(BICL, 1);axis([.5 Lmax+1.5 0 MM])set(gca, 'XTick', 1:Lmax+1)set(gca, 'XTickLabel', 0:Lmax)subplot(4, 2, 5);bar(AICM, 1);axis([.5 Mmax+0.5 0 MM])subplot(4, 2, 6);bar(AICL, 1);axis([.5 Lmax+1.5 0 MM])set(gca, 'XTick', 1:Lmax+1)set(gca, 'XTickLabel', 0:Lmax)subplot(4, 2, 7);bar(GICM, 1);axis([.5 Mmax+0.5 0 MM])subplot(4, 2, 8);bar(GICL, 1);axis([.5 Lmax+1.5 0 MM])set(gca, 'XTick', 1:Lmax+1)set(gca, 'XTickLabel', 0:Lmax)mean= [ sum((1:Mmax)'.*MDLM)/MM; sum((1:Mmax)'.*BICM)/MM; sum((1:Mmax)'.*AICM)/MM; sum((1:Mmax)'.*GICM)/MM; sum((1:Mmax)'.*MICM)/MM; sum((1:Mmax)'.*mdlM)/MM; sum((1:Mmax)'.*bicM)/MM; sum((1:Mmax)'.*aicM)/MM; sum((1:Mmax)'.*gicM)/MM; sum((1:Mmax)'.*micM)/MM; ];meanMn= [meanMn mean];mean= [ sum((0:Lmax).*MDLL)/MM; sum((0:Lmax).*BICL)/MM; sum((0:Lmax).*AICL)/MM; sum((0:Lmax).*GICL)/MM; sum((0:Lmax).*MICL)/MM; sum((0:Lmax).*mdlL)/MM; sum((0:Lmax).*bicL)/MM; sum((0:Lmax).*aicL)/MM; sum((0:Lmax).*gicL)/MM; sum((0:Lmax).*micL)/MM; ];meanLn= [meanLn mean];mse= [ 10*log10((sum(((1:Mmax)'-M).^2.*MDLM) + sum(((0:Lmax)-L).^2.*MDLL))/MM/(M^2+L^2)); 10*log10((sum(((1:Mmax)'-M).^2.*BICM) + sum(((0:Lmax)-L).^2.*BICL))/MM/(M^2+L^2)); 10*log10((sum(((1:Mmax)'-M).^2.*AICM) + sum(((0:Lmax)-L).^2.*AICL))/MM/(M^2+L^2)); 10*log10((sum(((1:Mmax)'-M).^2.*GICM) + sum(((0:Lmax)-L).^2.*GICL))/MM/(M^2+L^2)); 10*log10((sum(((1:Mmax)'-M).^2.*MICM) + sum(((0:Lmax)-L).^2.*MICL))/MM/(M^2+L^2)); 10*log10((sum(((1:Mmax)'-M).^2.*mdlM) + sum(((0:Lmax)-L).^2.*mdlL))/MM/(M^2+L^2)); 10*log10((sum(((1:Mmax)'-M).^2.*bicM) + sum(((0:Lmax)-L).^2.*bicL))/MM/(M^2+L^2)); 10*log10((sum(((1:Mmax)'-M).^2.*aicM) + sum(((0:Lmax)-L).^2.*aicL))/MM/(M^2+L^2)); 10*log10((sum(((1:Mmax)'-M).^2.*gicM) + sum(((0:Lmax)-L).^2.*gicL))/MM/(M^2+L^2)); 10*log10((sum(((1:Mmax)'-M).^2.*micM) + sum(((0:Lmax)-L).^2.*micL))/MM/(M^2+L^2)); ]; MSEn= [MSEn mse]; drawnowendN= 256;L= 2;MSEm= [];meanMm= [];meanLm= [];for M= 2:5 exp= 'ord'; resultname= sprintf('data/%04d/%s%d%d%s%s.mat', N, exp, M, L, re_im, mo_no) load(resultname) MDLM= sum(MDLhist.').';MDLL= sum(MDLhist); BICM= sum(BIChist.').';BICL= sum(BIChist); AICM= sum(AIChist.').';AICL= sum(AIChist); GICM= sum(GIChist.').';GICL= sum(GIChist); MICM= sum(MIChist.').';MICL= sum(MIChist); mdlM= sum(mdlhist.').';mdlL= sum(mdlhist); bicM= sum(bichist.').';bicL= sum(bichist); aicM= sum(aichist.').';aicL= sum(aichist); gicM= sum(gichist.').';gicL= sum(gichist); micM= sum(michist.').';micL= sum(michist);figure(M);clf;subplot(4, 2, 1);bar(MDLM, 1);axis([.5 Mmax+0.5 0 MM])subplot(4, 2, 2);bar(MDLL, 1);axis([.5 Lmax+1.5 0 MM])set(gca, 'XTick', 1:Lmax+1)set(gca, 'XTickLabel', 0:Lmax)subplot(4, 2, 3);bar(BICM, 1);axis([.5 Mmax+0.5 0 MM])subplot(4, 2, 4);bar(BICL, 1);axis([.5 Lmax+1.5 0 MM])set(gca, 'XTick', 1:Lmax+1)set(gca, 'XTickLabel', 0:Lmax)subplot(4, 2, 5);bar(AICM, 1);axis([.5 Mmax+0.5 0 MM])subplot(4, 2, 6);bar(AICL, 1);axis([.5 Lmax+1.5 0 MM])set(gca, 'XTick', 1:Lmax+1)set(gca, 'XTickLabel', 0:Lmax)subplot(4, 2, 7);bar(GICM, 1);axis([.5 Mmax+0.5 0 MM])subplot(4, 2, 8);bar(GICL, 1);axis([.5 Lmax+1.5 0 MM])?? 快捷鍵說明
復(fù)制代碼
Ctrl + C
搜索代碼
Ctrl + F
全屏模式
F11
切換主題
Ctrl + Shift + D
顯示快捷鍵
?
增大字號
Ctrl + =
減小字號
Ctrl + -