MATLAB實現(xiàn)LD算法進行AR估計
利用給定的一組樣本數(shù)據(jù)估計一個平穩(wěn)隨機信號的功率譜密度稱為功率譜估計，又稱譜估計。譜估計的方法可以分成經(jīng)典譜估計和現(xiàn)代譜估計。
經(jīng)典譜估計又稱為非參數(shù)化的譜估計，分為直接法和間接法。直接法是指直接計算樣本數(shù)據(jù)的傅里葉變換，即獲取頻譜，然后計算頻譜和其共軛的乘積，就得到功率譜；間接法是指先計算樣本數(shù)據(jù)的自相關(guān)函數(shù)，然后計算自相關(guān)函數(shù)的傅里葉變換，即得到功率譜。經(jīng)典譜估計存在很多的缺陷，主要原因是對數(shù)據(jù)加窗時默認在窗外未觀測到的數(shù)據(jù)的自相關(guān)系數(shù)為 0，這顯然是不切實際的；此外樣本數(shù)據(jù)是有限長的，而經(jīng)典譜估計往往需要較長的數(shù)據(jù)才能獲得較好性能，而且加窗函數(shù)也容易造成譜的模糊，因此我們需要參數(shù)化的譜估計，也就是現(xiàn)代譜估計。
現(xiàn)代譜估計中一類常用的模型是AR模型估計。
AR估計原理
AR 模型的表達式為：
上式兩邊同時作 Z 變換，得到：
AR 模型的原理就是根據(jù)觀測到的數(shù)據(jù)x(n)來估計 a i {a_i} ai?，從而估計輸入信號的功率譜。具體來說就是將輸入信號x(n)看成一個均值為 0，方差為 σ 2 \sigma^2 σ2的高斯噪聲通過信道得到的；先通過觀測數(shù)據(jù)估計信道參數(shù) a i {a_i} ai?，再結(jié)合 z 變換公式來計算輸入信號的功率譜。
AR 過程的線性預(yù)測和 Levinson-Durbin 算法
求解上述Yule-Walker方程的常用算法是 Levinson-Durbin 算法，具體推導(dǎo)過程如下：p-1 階的 Yule-Walker 方程如下：
代碼實現(xiàn)
算法步驟：
根據(jù) Levinson-Durbin 算法進行遞推，來計算 AR 模型的各階參數(shù)。采用最終預(yù)測誤差 FPE 法來確定 AR 模型階數(shù) p。
-
初始化：
-
迭代計算：
-
若k>p或者誤差小于閾值，返回各階次參數(shù)
-
對比估計信號和實際信號，觀察并分析結(jié)果
代碼：
LD_AR function函數(shù)存放在LD_AR.m文件中：文章來源：http://www.zghlxwxcb.cn/news/detail-768116.html
function [R,order, params, var] = LD_AR(x)
N = length(x)-1;
%% 自相關(guān)矩陣
R = zeros(N, 1);
for n = 1:N
R(n) = x(1:N-n+1) * (x(n:N))' /N;
end
%% LD 算法
a = zeros(N+1,N+1);
rho = zeros(1,N+1);
FPE = zeros(1,N+1);
rho(1) = R(1);
a(1,1) = -R(2)/R(1);
rho(2) = rho(1)*(1-abs(a(1,1))^2);
FPE(1) = (N+1)/(N-1) * rho(2);
for k=2:N-1
S = 0;
for m = 1:k-1
S = S + a(k-1,m) * R(k-m+1);
end
a(k,k) = -(R(k+1)+S)/rho(k);
for i = 1:k-1
a(k,i) = a(k-1,i) + a(k,k) * a(k-1,k-i);
end
rho(k+1) = rho(k)*(1-a(k,k)^2);
FPE(k) = (N+k)/(N-k) * rho(k+1);
end
%% 利用FPE確定AR模型階數(shù)
min = FPE(1);
disp 'min:'
disp(min)
for k = 2:N-1
if abs(FPE(k))<abs(min)
min = FPE(k);
order = k;
end
end
params = [1, a(order,1:order)];
var = abs(rho(order+1));
end
main.m文件：文章來源地址http://www.zghlxwxcb.cn/news/detail-768116.html
clear;
clc;
%% 產(chǎn)生一組樣本點
awgn_noise = randn(1,1000);
a_model = [1 -0.6 0.4 -0.3];
% a_model = [1 -0.4 0.2 -0.2 0.4];
x = filter(1, a_model, awgn_noise);
%% LD算法估計AR參數(shù)
[R,order, a, var] = LD_AR(x)
%% 原始信號 s
number=1000;
w = linspace(-pi,pi,number);
s = zeros([1,number]);
for m = 1:number
c = w(m);
kk = [1: length(a_model)-1];
s(m) = 1 /(abs(1+ a_model(2:end)* exp(-1i*c*kk')))^2;
end
%% 估計信號 s_hat
s_hat = zeros([1,number]);
for m = 1:number
c = w(m);
kk = [1: length(a)-1];
s_hat(m) = 1 /(abs(1+ a(2:end)* exp(-1i*c*kk')))^2;
end
%% 繪制頻譜
figure;
set(gcf,'position',[700, 200, 1400, 350])
subplot(1,3,1);
plot(w,s, 'b');
grid on
set(gca,'Xtick',[-pi,-pi/2,0, pi/2,pi]);
set(gca,'XTickLabel',{'-pi' '-pi/2' '0' 'pi/2' 'pi'});
title('原始信號功率譜','fontsize', 12);
xlabel('頻率','fontsize', 12);
ylabel('功率','fontsize', 12);
ylim([0,7]);
subplot(1,3,2);
plot(w,s_hat, 'r-');
grid on
set(gca,'Xtick',[-pi,-pi/2,0, pi/2,pi]);
set(gca,'XTickLabel',{'-pi' '-pi/2' '0' 'pi/2' 'pi'});
title('LD算法進行AR估計功率譜', 'fontsize', 12);
xlabel('頻率','fontsize', 12);
ylabel('功率','fontsize', 12);
ylim([0,7]);
subplot(1,3,3);
plot(w,s,'b');
hold on
plot(w,s_hat, 'r--');
hold off
grid on
set(gca,'Xtick',[-pi,-pi/2,0, pi/2,pi]);
set(gca,'XTickLabel',{'-pi' '-pi/2' '0' 'pi/2' 'pi'});
title('原始信號功率譜和LD算法進行AR估計比較', 'fontsize', 12);
xlabel('頻率','fontsize', 12);
ylabel('功率','fontsize', 12);
legend('原始信號','AR估計信號')
ylim([0,7]);
% clc;
% clear;
%
% % 構(gòu)造加入高斯白噪聲的樣本點函數(shù)
% f1=0.4;
% f2=0.2;
% N=input('產(chǎn)生的樣本點個數(shù)N:');
% n=1:N;
% wn=randn(1,N);% 產(chǎn)生高斯白噪聲
% xn=sin(2*pi*f1*n+pi/3)+10*sin(2*pi*f2*n+pi/4);
% yn=xn+wn;
%
% Rx=xcorr(yn,'biased');
%
% k=1:N;
% f=k/N;
% Sx=abs(fft(Rx,N));
% plot(f,Sx);
% title('功率譜');
% xlabel('f');
% ylabel('Sx');
%
% % 根據(jù)自相關(guān)函數(shù)估計f1,f2
% sxk=zeros(1,N/2);
% for k=1:N/2
% sxk(k)=Sx(k);
% end
%
% [pks,locs]=findpeaks(sxk,'SortStr','descend');
% locs(1)
% locs(2)
% for i=1:2
% fprintf('預(yù)測結(jié)果：f%d=%f\n',i,locs(i)/N);
% end
到了這里，關(guān)于現(xiàn)代信號處理實驗：MATLAB實現(xiàn)LD算法進行AR估計的文章就介紹完了。如果您還想了解更多內(nèi)容，請在右上角搜索TOY模板網(wǎng)以前的文章或繼續(xù)瀏覽下面的相關(guān)文章，希望大家以后多多支持TOY模板網(wǎng)！
