levension.m

来自「列文森-杜宾算法」· M 代码 · 共 110 行

s1=wavread('8a.wav');
s=s1(3800:4055);
d=hanning(256);
w=s.*d;
t=fft(w,256);
h=fftshift(t);
nn=log(abs(h));
gg=nn(128:255);
p=4;
for j=1:p
    r(j)=w(j:256)'*w(1:256-j+1);
end
[A E K]=levinson(r,p);
a1=fft(A,256);
hh=log(p*E)-log(abs(a1));
hh1=hh(1:128);
p=8;
for j=1:p
    r(j)=w(j:256)'*w(1:256-j+1);
end
[A E K]=levinson(r,p);
a1=fft(A,256);
hh=log(p*E)-log(abs(a1));
hh2=hh(1:128);
f=[0:127]*8000/256;

p=12;
for j=1:p
    r(j)=w(j:256)'*w(1:256-j+1);
end
[A E K]=levinson(r,p);
a1=fft(A,256);
hh=log(p*E)-log(abs(a1));
hh3=hh(1:128);
p=16;
for j=1:p
    r(j)=w(j:256)'*w(1:256-j+1);
end
[A E K]=levinson(r,p);
a1=fft(A,256);
hh=log(p*E)-log(abs(a1));
hh4=hh(1:128);
p=20;
for j=1:p
    r(j)=w(j:256)'*w(1:256-j+1);
end
[A E K]=levinson(r,p);
a1=fft(A,256);
hh=log(p*E)-log(abs(a1));

hh5=hh(1:128);

p=30;
for j=1:p
    r(j)=w(j:256)'*w(1:256-j+1);
end
[A E K]=levinson(r,p);
a1=fft(A,256);
hh=log(p*E)-log(abs(a1));
hh6=hh(1:128);

subplot(3,3,1)
plot(f,gg);
hold on
plot(f,hh1)
xlabel('Frequency/Hz');
ylabel('Magnitude/dB');
text(3000,3,'p=4');
subplot(3,3,2)
plot(f,gg);
hold on
plot(f,hh2)
xlabel('Frequency/Hz');
ylabel('Magnitude/dB');
text(2500,3,'p=8');
subplot(3,3,3)
plot(f,gg);
hold on
plot(f,hh3)
xlabel('Frequency/Hz');
ylabel('Magnitude/dB');
text(2500,3,'p=12');
subplot(3,3,4)
plot(f,gg);
hold on
plot(f,hh4)
xlabel('Frequency/Hz');
ylabel('Magnitude/dB');
text(2500,3,'p=16');
subplot(3,3,5)
plot(f,gg);
hold on
plot(f,hh5)
xlabel('Frequency/Hz');
ylabel('Magnitude/dB');
text(2500,3,'p=20');
subplot(3,3,6)
plot(f,gg);
hold on
plot(f,hh6)
xlabel('Frequency/Hz');
ylabel('Magnitude/dB');
text(2500,3,'p=30');