1. AR模型概念觀
AR模型是一種線性預測，即已知N個數據，可由模型推出第N點前面或后面的數據（設推出P點），所以其本質類似於插值，其目的都是為了增加有效數據，只是AR模型是由N點遞推，而插值是由兩點（或少數幾點）去推導多點，所以AR模型要比插值方法效果更好。
數字信號處理功率譜估計方法分經典功率譜估計和現代功率譜估計，現代功率譜估計以參數模型功率譜估計為代表，參數功率譜模型如下：
u(n) ——> H(z) ——> x(n)
參數模型的基本思路是：
—— 參數模型假設研究過程是由一個輸入序列u(n)激勵一個線性系統H(z)的輸出。
—— 由假設參數模型的輸出x(n)或其自相關函數
來估計H(z)的參數
—— 由H(z)的參數估計x(n)的功率譜
因此，參數模型功率譜的求解有兩步：
（1）H(z)模型參數估計
（2）依據模型參數求功率譜
AR模型（自回歸模型，Auto Regression Model）是典型的現代參數功模型。其定義為
其中，輸入設定為方差為
的白噪聲序列，ak是模型的參數，p是模型的階數，Px為x(n)功率譜，也即本文要求解的目標。
AR模型是一個全極點模型，“自回歸”的含義是：現在的輸出是現在的輸入和過去p個輸出的加權和。
現在我們希望建立AR參數模型和x(n)的自相關函數的關系，也即AR模型的正則方程：
上面的正則方程也稱Yule-Walker方程，其中的rx為自相關函數。由方程可以看出，一個p階的AR模型有p+1個參數（
）。
通過推導可以發現，AR模型與線性預測器是等價的，AR模型是在最小平方意義上對數據的擬合。
2. AR模型參數求解——Levinson-Durbin Algorithm
定義
為p階AR模型在m階次時的第k個系數，k=1,2,...,m。定義
為m階系統時的，這也是線性預測器中前向預測的最小誤差功率。此時，一階AR模型時有
我們定義初始時
，則
由PART1中矩陣的對稱性質，將上面的公式推廣到高階AR模型，可以推導出Levinson-Durbin遞推算法：
Levinson-Durbin遞推算法從低階開始遞推，，給出了每一階次時所有參數，。這一特點有利於我們選擇合適的AR模型階次。
因為
必須大於0，由式
知
，如果
，遞推應該停止。
到此，選擇最佳階次的參數代入到
中，求得功率譜。
3. matlab實現
matlab工具箱中提供了現成的函數實現AR模型功率譜計算。參考[2]，我們將內容摘錄如下：
AR模型的譜估計是現代譜估計的主要內容。
1．AR模型的Yule—Walker方程和Levinson-Durbin遞推算法：在MATLAB中，函數levinson和aryule都采用Levinson-Durbin遞推算法來求解AR模型的參數a1,a2,……,ap及白噪聲序列的方差，只是兩者的輸入參數不同，它們的格式為：
A=LEVINSON(R,ORDER) A=ARYULE(x,ORDER)
兩函數均為定階ORDER的求解，但是函數levinson的輸入參數要求是序列的自相關函數，而函數aryule的輸入參數為采樣序列。
下面語句說明函數levinson和函數aryule的功能是相同的：
例子：
randn('seed',0)
a=[1 0.1 0.2 0.3 0.4 0.5];
x=impz(1,a,20)+randn(20,1)/20;
r=xcorr(x,'biased');
r(1:length(x)-1)=[];
A=levinson(r,5)
B=aryule(x,5)
2．Burg算法：
格式為：A=ARBURG(x,ORDER); 其中x為有限長序列，參數ORDER用於指定AR模型的階數。以上面的例子為例：
randn('seed',0)
a=[1 0.1 0.2 0.3 0.4 0.5];
x=impz(1,a,20)+randn(20,1)/20;
A=arburg(x,5)
3.改進的協方差法：
格式為：A=ARMCOV(x,ORDER); 該函數用來計算有限長序列x(n)的ORDER階AR模型的參數。例如：輸入下面語句：
randn('seed',0)
a=[1 0.1 0.2 0.3 0.4 0.5];
x=impz(1,a,20)+randn(20,1)/20;
A=armcov(x,5)
AR模型階數P的選擇：
AR模型階數P一般事先是不知道的，需要事先選定一個較大的值，在遞推的過程中確定。在使用Levinson—Durbin遞推方法時，可以給出由低階到高階的每一組參數，且模型的最小預測誤差功率Pmin（相當於白噪聲序列的方差）是遞減的。直觀上講，當預測誤差功率P達到指定的希望值時，或是不再發生變化時，這時的階數即是應選的正確階數。
因為預測誤差功率P是單調下降的，因此，該值降到多少才合適，往往不好選擇。比較常見的准則是：
最終預測誤差准則：FPE(r)=Pr{[N+(r+1)]/ [N-(r+1)]}
信息論准則：AIC(r)=N*log(Pr)+2*r
上面的N為有限長序列x(n)的長度，當階數r由1增加時，FPE(r) 和AIC(r)都將在某一r處取得極小值。將此時的r定為最合適的階數p。
MATLAB中AR模型的譜估計的函數說明：
1． Pyulear函數：
功能：利用Yule--Walker方法進行功率譜估計.
格式： Pxx=Pyulear(x,ORDER,NFFT)
[Pxx,W]=Pyulear(x,ORDER,NFFT)
[Pxx,W]=Pyulear(x,ORDER,NFFT,Fs)
Pyulear(x,ORDER,NFFT,Fs,RANGE,MAGUNITS)
說明：Pxx =Pyulear(x,ORDER,NFFT)中，采用Yule--Walker方法估計序列x的功率譜，參數ORDER用來指定AR模型的階數，NFFT為FFT算法的長度，默認值為256，若NFFT為偶數，則Pxx為（NFFT/2 + 1）維的列矢量，若NFFT為奇數，則Pxx為（NFFT + 1）/2維的列矢量；當x為復數時，Pxx長度為NFFT。
[Pxx,W]=Pyulear(x,ORDER,NFFT)中，返回一個頻率向量W.
[Pxx,W]=Pyulear(x,ORDER,NFFT,Fs)中，可以在F向量得到功率譜估計的頻率點，Fs指定采樣頻率。
Pyulear(x,ORDER,NFFT,Fs,RANGE,MAGUNITS)中，直接畫出功率譜估計的曲線圖。
2． Pburg函數：
功能：利用Burg方法進行功率譜估計。
格式：Pxx=Pburg(x,ORDER,NFFT)
[Pxx,W]=Pburg(x,ORDER,NFFT)
[Pxx,W]=Pburg(x,ORDER,NFFT,Fs)
Pburg(x,ORDER,NFFT,Fs,RANGE,MAGUNITS)
說明：Pburg函數與Pyulear函數格式相同，只是計算AR模型時所采用的方法不同，因此格式可以參照Pyulear函數。
3． Pcov函數：
功能：利用協方差方法進行功率譜估計。
格式：Pxx=Pcov(x,ORDER,NFFT)
[Pxx,W]=Pcov(x,ORDER,NFFT)
[Pxx,W]=Pcov(x,ORDER,NFFT,Fs)
Pcov(x,ORDER,NFFT,Fs,RANGE,MAGUNITS)
說明：Pcov函數采用協方差法估計AR模型的參數，然后計算序列x的功率譜。協方差法與改進的協方差法相比，前者僅令前向預測誤差為最小，其他步驟是一樣的。：Pcov函數與Pyulear函數格式相同，只是計算AR模型時所采用的方法不同，因此格式可以參照Pyulear函數.
4.Pmcov:
功能：利用改進的協方差方法進行功率譜估計。
格式：Pxx=Pmcov(x,ORDER,NFFT)
[Pxx,W]=Pmcov(x,ORDER,NFFT)
[Pxx,W]=Pmcov(x,ORDER,NFFT,Fs)
Pmcov(x,ORDER,NFFT,Fs,RANGE,MAGUNITS)
例如：輸入下面語句：
figure 8.10--8.11
Fs=1000; %采樣頻率
n=0:1/Fs:3;
xn=cos(2*pi*n*200)+randn(size(n));
%設置參數
order=20;
nfft=1024;
%Yule-Walker方法
figure(1)
pyulear(xn,order,nfft,Fs);
%Burg方法
figure(2)
pburg(xn,order,nfft,Fs);
%協方差法
figure(3)
pcov(xn,order,nfft,Fs);
%改進協方差方法
figure(4)
pmcov(xn,order,nfft,Fs);
AR譜的分辨率：
經典譜估計的分辨率反比與信號的有效長度,但是現代譜估計的分辨率可以不受此限制. 這是因為對於給定的N點有限長序列x(n)，雖然其估計出的相關函數也是有限長的，但是現代譜估計的一些方法隱含着數據和自相關函數的外推，使其可能的長度超過給定的長度，因而AR譜的分辨率較高。
例如：序列x(n)由兩個正鉉信號組成，其頻率分別為f1=20Hz和f2=21Hz,並含有一定的噪聲量。試分別用周期圖法，Burg方法與改進的協方差法估計信號的功率譜，且AR模型的階數取30和50兩種情況討論。
上面的例子可以通過下面程序實現：
Fs=200;
n=0:1/Fs:1;
xn=sin(2*pi*20*n)+sin(2*pi*21*n)+0.1*randn(size(n));
window=boxcar(length(xn));
nfft=512;
[Pxx,f]=periodogram(xn,window,nfft,Fs);
figure(1)
plot(f,10*log10(Pxx)),grid
xlabel('Frequency(Hz)')
ylabel('Power Spectral Density(dB/Hz)')
title('Periodogram PSD Estimate')
order1=30;
order2=50;
figure(2)
pburg(xn,order1,nfft,Fs)
figure(3)
pburg(xn,order2,nfft,Fs)
figure(4)
pmcov(xn,order1,nfft,Fs)
figure(5)
pmcov(xn,order1,nfft)
LMS自適應濾波器是使濾波器的輸出信號與期望響應之間的誤差的均方值為最小，因此稱為最小均方（LMS）自適應濾波器。其原理及推導見http://download.csdn.net/source/432206
function [yn,W,en]=LMS(xn,dn,M,mu,itr) % LMS(Least Mean Squre)算法 % 輸入參數: % xn 輸入的信號序列 (列向量) % dn 所期望的響應序列 (列向量) % M 濾波器的階數 (標量) % mu 收斂因子(步長) (標量) 要求大於0,小於xn的相關矩陣最大特征值的倒數 % itr 迭代次數 (標量) 默認為xn的長度,M<itr<length(xn) % 輸出參數: % W 濾波器的權值矩陣 (矩陣) % 大小為M x itr, % en 誤差序列(itr x 1) (列向量) % yn 實際輸出序列 (列向量) % 參數個數必須為4個或5個 if nargin == 4 % 4個時遞歸迭代的次數為xn的長度 itr = length(xn); elseif nargin == 5 % 5個時滿足M<itr<length(xn) if itr>length(xn) | itr<M error('迭代次數過大或過小!'); end else error('請檢查輸入參數的個數!'); end % 初始化參數 en = zeros(itr,1); % 誤差序列,en(k)表示第k次迭代時預期輸出與實際輸入的誤差 W = zeros(M,itr); % 每一行代表一個加權參量,每一列代表-次迭代,初始為0 % 迭代計算 for k = M:itr % 第k次迭代 x = xn(k:-1:k-M+1); % 濾波器M個抽頭的輸入 y = W(:,k-1).' * x; % 濾波器的輸出 en(k) = dn(k) - y ; % 第k次迭代的誤差 % 濾波器權值計算的迭代式 W(:,k) = W(:,k-1) + 2*mu*en(k)*x; end % 求最優時濾波器的輸出序列 yn = inf * ones(size(xn)); for k = M:length(xn) x = xn(k:-1:k-M+1); yn(k) = W(:,end).'* x; end
LMS函數的一個實例：
%function main() close all % 周期信號的產生 t=0:99; xs=10*sin(0.5*t); figure; subplot(2,1,1); plot(t,xs);grid; ylabel('幅值'); title('it{輸入周期性信號}'); % 噪聲信號的產生 randn('state',sum(100*clock)); xn=randn(1,100); subplot(2,1,2); plot(t,xn);grid; ylabel('幅值'); xlabel('時間'); title('it{隨機噪聲信號}'); % 信號濾波 xn = xs+xn; xn = xn.' ; % 輸入信號序列 dn = xs.' ; % 預期結果序列 M = 20 ; % 濾波器的階數 rho_max = max(eig(xn*xn.')); % 輸入信號相關矩陣的最大特征值 mu = rand()*(1/rho_max) ; % 收斂因子 0 < mu < 1/rho [yn,W,en] = LMS(xn,dn,M,mu); % 繪制濾波器輸入信號 figure; subplot(2,1,1); plot(t,xn);grid; ylabel('幅值'); xlabel('時間'); title('it{濾波器輸入信號}'); % 繪制自適應濾波器輸出信號 subplot(2,1,2); plot(t,yn);grid; ylabel('幅值'); xlabel('時間'); title('it{自適應濾波器輸出信號}'); % 繪制自適應濾波器輸出信號,預期輸出信號和兩者的誤差 figure plot(t,yn,'b',t,dn,'g',t,dn-yn,'r');grid; legend('自適應濾波器輸出','預期輸出','誤差'); ylabel('幅值'); xlabel('時間'); title('it{自適應濾波器}');