%----file read
%soundsc(sp.in,8000,16)
soundscってやつでsp.inってやつの音が聞けます
ちなみにsp.inのように変数に . を使うのってよくないそうなので避けた方がいいです。
%B =fir1(500,1/4);
%sp.in= conv(sp.in,B);
ここでローパスフィルタをかけるつもりでした。
経験的にfs>8khzではあまり意味がないようです。
こっからあとは授業のソースをほぼ丸写し。
ds = filter([1,-0.9],1,sp.in); % pre-emphasis (1st order high-pass filtering)
period = zeros( 1 , nseg );
inmax=max(sp.in);
figure(100);
subplot(5,1,1)
smax = max([max(sp.in),-min(sp.in)]);
plot(1:len,sp.in(is:ie))
axis([1 len -smax smax]);
ylabel('Speech signal')
% segmentation and LPC analysis of the "differentiated" speech
a = zeros(1,order+1);
ここにはフィルタ係数が入ります。
Cを使う人、callocみたいなものをやってるとおもってください。
res = zeros(1,length(sp.in)); % residual signal will be stored here
これは残差信号です。
長さが窓の長さ(Nf)と一致していることを覚えておいてください。
w = hamming(Nf); % hamming window of length N
for l=1:nseg
ここからループをまわします。
これは音声波形をある感覚毎に区切っていくことを意味します。
区切られたフレームを分析フレームと以降いいます。
iss = is+Ns*(l-1);
分析フレームのはじめのインデックスをissと置いています
x = ds(iss:iss+Nf-1)'.*w'; % segmentation using the window
a1 = a; % saving the coefficients of the previous frame
ここで前フレームのフィルタ係数を保存(補間じゃない)してます。
a = lpc(x,order); % LPC analysis
[h,wf] = freqz(1,a,nfft,'whole'); % frequency response of the filter
spec = abs(h); % magnitude spectrum
smax = max(spec);
spec = 20*log10(spec/smax); % log magnitude spectrum
subplot(5,1,2)
spec = (Nf/2+Ns*(l-1))+Ns+Ns*spec/16; % 16 is a magic number
ところでmagic numberってなんなのでしょう?だれか教えて。
f = linspace(0,1,nfft/2);
plot(spec(1:nfft/2),f); hold on;
axis([1 len 0 1]);
ylabel('Normalized frequency');
% inverse filtering of the "original" speech
ここからフィルタ係数の補間に入ります
if (l==1) %%% 1st analysis frame %%%
aa(1:order+1) = a(order+1:-1:1); % inversion of the vector order
for m=1:Nf/2
res(m) = dot(aa,sp.in(is+m-order-1:is+m-1)); % inverse filtering
end
else %%% 2nd and the succeedin asnalysis frames %%%
k1 = poly2rc(a1); % conversion of the spectral paramaters
k = poly2rc(a);
for m=1:Ns % processing for the individual sample
kk = k1+(k-k1).*(m-1)/(Ns-1); % spectral interpolation
aa = rc2poly(kk); % conversion of the spectral paramaters
aa(1:order+1) = aa(order+1:-1:1); % inversion of the vector order
ix = Nf/2+Ns*(l-2)+m;
res(ix) = dot(aa,sp.in(is+ix-order-1:is+ix-1)); %inverse filetering
end
end
ここまでの流れを各フレームに対しておこないます。
end
ここで注意したいのは
subplot(5,1,2)
でプロットしたものは補間する前のフィルタ係数なんです。
というのもうちの班ではあまり必要なデータではなかったのでそうしています。
補間後のフィルタ係数が出力したければ、
「ここまでの流れを各フレームに対しておこないます。」
と記述した辺りに
[h,wf] = freqz(1,a,nfft,'whole'); % frequency response of the filter
spec = abs(h); % magnitude spectrum
smax = max(spec);
spec = 20*log10(spec/smax); % log magnitude spectrum
subplot(5,1,2)
spec = (Nf/2+Ns*(l-1))+Ns+Ns*spec/16; % 16 is a magic number
のaをaaにしたらいいはず。
subplot(5,1,3)
%ls = N+L*(nseg-1);
res = real(res);
rmax = (max([max(res),-min(res)]));
plot(1:len,res(1:len))
axis([1 len -rmax rmax]);
ylabel('Residual signal');
で残差信号のプロット。
ここでちゃんと波形がでているか確認してください。
ここからは別にしなくてもいいんですけど、
残差信号をLPCかけて、残差信号のスペクトル抱絡を出力しようとしています。
残差信号にLPCって感覚がちょっとなれないですよね(´ⅴ`;
% segmentation and LPC analysis of the residual signal
for l=1:nseg
iss = 1+Ns*(l-1);
x = res(iss:iss+Nf-1).*w'; % segmentation using the window
a = lpc(x,order); % LPC analysis
[h,wf] = freqz(1,a,nfft,'whole'); % frequency response of the filter
spec = abs(h); % magnitude spectrum
smax = max(spec);
spec = 20*log10(spec/smax); %log magnitude spectrum
subplot(5,1,4)
spec = (Nf/2+Ns*(l-1))+Ns+Ns*spec/16; % 16 is a magic number
f = linspace(0,1,nfft/2);
plot(spec(1:nfft/2),f); hold on;
axis([1 len 0 1]);
ylabel('Normalized frequency');
end