功率譜估計 —— AR(Autoregressive)

function [A,Sgm2] = AR(x,P)% x - 輸入數據，列向量% P - AR 階數,標量% A - AR 係數[a1 a2 ...aP]，行向量;A = zeros(1,P); % 初始化參數矩陣，行向量Rx = Rxx(x,P); % 無偏自相關估計，size = P+1% p = 1 時候的Yule-Walker方程解,注意，MATLAB下標從1開始A(1) = - Rx(1+1)/Rx(0+1);Sgm2 = Rx(0+1)*(1 - A(1)^2);% p = 1->P 時候的遞推,注意，MATLAB下標從1開始for p = 1:P-1 k = 1:p; K = -(Rx((p+1)+1) + A(k)*Rx((p+1-k)+1))/Sgm2; Sgm2 = Sgm2 * (1 - K*K); A(k) = A(k) + K * A(p+1-k); A(p+1) = K;end

所用函數二：

function Rx = Rxx(x,P)% x - 自變數% P - 自相關數目% Rxx[m] = sum(x[n]*x[n+m]),n = 0->inf,m = 0->P;N = length(x);sizex = size(x);if sizex(1)==1 Rx = zeros(1,P+1); % x 為行向量那麼Rx也為行向量 for m = 0:P Rx(m+1) = x(m+1:N)*x(1:N-m); end elseif sizex(2)==1 Rx = zeros(P+1,1); % x 為列向量那麼Rx也為列向量 for m = 0:P Rx(m+1) = x(m+1:N)*x(1:N-m); end end

所用函數三：

function [H,f] = ar2fhz(A,fs,df)% A - AR 的係數[a1 a2 ... ap]，行向量;% fs - 指定的頻率範圍% dim - 解析度,角度；% H - 功率譜密度函數，行向量% df - 頻率解析度，標量P = length(A);% Hz解析度->弧度解析度dw = df/fs*2*pi; w = 0:dw:pi;z = exp(-j*w);H = zeros(size(w));for i = 1:P H = H + A(i)*z.^i;endf = w/2/pi*fs;H = 1./(1 + H);H = abs(H).^2;

測試

ticfs = 600; % 採樣率t = 0:1/fs:1;f1 = 100; % 頻率100赫茲f2 = 160; % 頻率160赫茲f3 = 230; % 頻率230赫茲x1 = sin(2*pi*f1*t); x2 = 1.2*sin(2*pi*f2*t); x3 = 1.5*sin(2*pi*f3*t); x = x1 + x2 + x3;df = 1; % 頻域間隔[a,~] = AR(x,40);[H,f] = ar2fhz(a,fs,df);plot(f,log(H),r)toc

結果

還行，還能用。。。
推薦閱讀：