由 LPC 參數到 LPCC 參數

01-28

&&【主體函數】&&

1. 由 LPC 參數轉換到 LPCC 參數函數

function h = LPC2LPCC(a,M)nn% 由AR模型導得LPC係數轉換為LPCC係數n% 1 n% H(z) = --------------------n% 1 - sum(a[k]*z^(-k))n% a - [a1 a2 a3 ... aP];n% M - 保留的 h(n) 係數個數nna = -a; % 在AR2LPC里分母項是+nP = length(a); % AR模型的階數nh = zeros(1,M); % lpcc的個數nh(1) = a(1); % 第一個n% 方法一nfor m = 2:Mn if m <= Pn k = 1:m-1;n h(m) = a(m) + sum((1-k/m).*a(k).*h(m-k));n elseif m > Pn k = 1:P;n h(m) = sum((1 - k/m).*a(k).*h(m-k)); n endnendnn% 方法二n% for m = 2:Mn% if m <= Pn% k = 1:m-1;n% h(m) = a(m) + sum((k/m).*a(m-k).*h(k));n% elseif m > Pn% k = m-P:m-1;n% h(m) = sum((k/m).*a(m-k).*h(k)); n% endn% endn

2. 由 Levinson - Durbin 遞推演算法得到 Yule - Walker 方程的解

function A = AR2LPC(x,P)nn% 根據AR模型輸出線型預測係數n% 1 n% H(z) = --------------------n% 1 + sum(a[k]*z^(-k))n% function A = ar2lpc(x,P)n% x - 輸入數據，列向量n% P - AR 階數,標量n% A - AR 係數[a1 a2 ...aP]，行向量;n% Sgm2 - 方差nnA = zeros(1,P); % 初始化參數矩陣，行向量nRx = Rxx(x,P); % 無偏自相關估計，size = P+1nn% p = 1 時候的Yule-Walker方程解,注意，MATLAB下標從1開始nA(1) = - Rx(1+1)/Rx(0+1);nSgm2 = Rx(0+1)*(1 - A(1)^2);nn% p = 1->P 時候的遞推,注意，MATLAB下標從1開始nfor p = 1:P-1n k = 1:p;n K = -(Rx((p+1)+1) + A(k)*Rx((p+1-k)+1))/Sgm2;n Sgm2 = Sgm2 * (1 - K*K);n n A(k) = A(k) + K * A(p+1-k);n A(p+1) = K;nendn

3. 自相關函數

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

#【用於驗證的函數】#

1. 由 LPC參數轉換到頻域

function [H,f] = LPC2HZ(A,Sgm2,fs,df)n% A - AR 的係數[a1 a2 ... ap]，行向量;n% Sgm2 - 方差n% fs - 指定的頻率範圍n% dim - 解析度,角度；n% H - 功率譜密度函數，行向量n% df - 頻率解析度，標量nnP = length(A);nndw = df/fs*2*pi; % Hz解析度->弧度解析度nw = 0:dw:pi;nz = exp(-j*w);nH = zeros(size(w));nfor i = 1:Pn H = H + A(i)*z.^i;nendnf = w/2/pi*fs;nH = 1./(1 + H);nH = Sgm2*abs(H).^2;n

2. 由單位脈衝響應轉化到頻域

function H = pulse2hz(h,fs,df)nn% H - 頻響n% h - 單位脈衝響應函數n% fs - 指定的頻率範圍n% df - 頻率解析度，標量nn% Hz解析度->弧度解析度ndw = df/fs*2*pi; nw = 0:dw:pi;nn = length(h);nz = exp(-j*w); nH = zeros(size(w));nfor i = 0:n-1n H = H + h(i+1).*z.^i;nendn

測試腳本 & 結果：

fs = 1023;nt = 0:1/1023:1;nx = sin(2*pi*100*t)+sin(2*pi*300*t);nN = length(x);nna = AR2LPC(x,10);n[H,F] = LPC2HZ(a,1,fs,1);nsubplot(211)nplot(F,log(H),b)naxis tightnnh = LPC2LPCC(a,100);nH = pulse2hz(h,fs,1);nsubplot(212)nplot(abs(H),r)naxis tightn

用兩個單頻信號模擬就這樣了。。。