人気ブログランキング | 話題のタグを見る

room415


SIFT アルゴ2 フィルタ係数の補間まで

%----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
by y0nda | 2005-07-20 20:41

<< SIFT アルゴ3 自己相関お...      SIFT アルゴリズム1 >>

日々取り組んでいることを書き留めています。
by y0nda
S M T W T F S
1 2
3 4 5 6 7 8 9
10 11 12 13 14 15 16
17 18 19 20 21 22 23
24 25 26 27 28 29 30
31
カテゴリ
お気に入りブログ
リンク
その他のジャンル
ファン
記事ランキング
ブログジャンル
画像一覧