%****************************************************************
%six?method?to?compute?spectrum
%?1:?classic?method
%?2::AR?model?method
%?3:?Pisarenko?method
%?4:?MUSIC?method
%?5:?ESPRIT?method
%?6:?TLS-ESPRIT?method
%?by?jianjian?2011/11/24
%****************************************************************
clear?all
clc
%sample?points?number
N?=?128;
%?runtime?number
runNum?=?20;
%?frequecy
f1?=?0.05;
f2?=?0.30;
f3?=?0.31;

n?=?[0:N-1]‘;


noiseVar?=?0.16;
%?mean?and?var?changed?with?noise
%?verify?how?many?complex?signals
%*********************
noiseAmp?=?sqrt（noiseVar）;
inputData?=?2*cos（2*pi*f1*n）?+?3*cos（2*pi*f2*n）?+?1.2*cos（2*pi*f3*n）?+?noiseAmp*randn（N1）;
M?=?N-100;?%?define?M
rx?=?corrmtx（inputDataM-1）;?%?corr?matrix
Rxx?=?rx‘*rx;
Lamda?=?abs（eig（Rxx））;
Lamda?=?sort（Lamda‘descend‘）;
mdl?=?zeros（M1）;?%?choose?the?smallest?mdl?to?verify?order
for?pTestCount?=?0:M-1
????A?=?（sum（Lamda（pTestCount+1:M1））/（M-pTestCount））^（M-pTestCount）;
????G?=?prod（Lamda（pTestCount+1:M1））;
????mdl（pTestCount+11）?=?-N*log10（G/A）+0.5*pTestCount*（2*M-pTestCount）*log10（N）;
end
err?=?0.03;?%?calculate?error?permitted
temp_mdl?=?diff（mdl/max（mdl））;
for?pTestCount?=?1:M
????if（abs（temp_mdl（pTestCount1））????????break;
????end
end
M?=?pTestCount;
pDimension?=?M-1;
%******************************

%?output?frequece
outputFreq?=?zeros（NrunNum5）;
FreqEval?=?zeros（pDimensionrunNum6）;
Feather?=?zeros（pDimension26）;


%?calculate
%***************************************
for?count?=?1:20
????%generate?noise?and?signal
????noise?=?noiseAmp*randn（N1）;
????inputData?=?2*cos（2*pi*f1*n）?+?3*cos（2*pi*f2*n）?+?1.2*cos（2*pi*f3*n）?+?noise;
????%?classic?method
????outputFreq（:count1）?=?abs（fft（inputData））.^2/N;
????????%?spectrum?evaluate
????temp_Spectrum?=?diff（outputFreq（:count1））;?%?diff?to?find?max?easily
????temp_Spectrum?=?temp_Spectrum/max（temp_Spectrum）;
????temp_Freq?=?zeros（pDimension1）;?%?iteration?to?find?spectrum?peak
????err?=?0.03;
????FCount?=?0;
????while（FCount?~=?pDimension）
????????FCount?=?0;
????????for?pTestCount?=?1:size（temp_Spectrum1）-1
????????????if?（temp_Spectrum（pTestCount1）>err?&&?temp_Spectrum（pTestCount+11）<-1*err）
????????????????FCount?=?FCount?+?1;
????????????????if（FCount?<=?pDimension）
????????????????????temp_Freq（FCount1）?=?pTestCount;
????????????????end
????????????end
????????end
????????if（FCount?>?pDimension）
????????????err?=?err?+?0.03;
????????else
????????????err?=err?-?0.03;
????????end
????????if（FCount?==?pDimension）
????????????break;
????????end
????end
????FreqEval（:count1）?=?temp_Freq/N;
????
????%?AR?method
????arObj?=?iddata（inputData）;
????p_test?=?zeros（N/101）;
????for?pTestCount?=?1:N/10???%calculate?AIC
????????ARPolyModel?=?ar（arObjpTestCount‘yw‘‘ppw‘）;
????????AIC?=?aic（ARPolyModel）;
????????p_test（pTestCount1）?=??AIC;
????end
????[AIC?pOrder]?=?min（p_test）;??%choose?the?smallest?AI