matlab版本频率估计

%****************************************************************
%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