目次
音のスペクトル解析を行った所では、振幅スペクトルと位相スペクトルについて
紹介して来ました。横軸が周波数で、縦軸が振幅や位相の成分でした。音のデータ
を特に音声データに限ると特徴的なスペクトルが現れます。特徴的に現れるスペク
トルの「音声のパラメータ」について今回は学びます。
6-2LPC(線形予測分析)
LPCについては火曜日のゼミの本で詳しくやりますので、今回は説明を省略します。
x(n)=g(n)*v(n)
この信号x(n)に離散フーリエ変換の絶対値の対数演算操作Dを掛けると、対数は掛け算を和の形に
するので以下の様になる。
D{x(n)}=D{g(n)*v(n)}=D{g(n)}+D{v(n)}
(注:)積が和に置き換えられるのは対数演算のためである。
このときD{x(n)}の逆離散フーリエ変換をケプストラムc(n)と呼ぶ。対数演算によりスペクトルの
高周波と低周波が分離され、逆離散フーリエ変換を施すと時間軸の低次のところにスペクトル低周波
の共振情報が、時間軸の高次のところに音源情報が現れる。(逆離散フーリエ変換を行うので、ケプス
トラム変数は時間の次元を持つ。)
次のMATLABプログラムを実行しなさい。またこのプログラムについて以下の質問に答えなさい。
============================================================================
フレーム表示のMOVIEを作るプログラム
============================================================================
課題が出来たら、答えを山本(eli@sys.wakayama-u.ac.jp)までメイルで送って下さい。
締め切りは12/7(木)17:00までです。
close all; % 画面クリア clear all; % データクリア [x,fs]=aiffread('kousyouAD16.aiff'); %データ読み込み t=[0:1/fs:length(x)/fs-1/fs]; %時間軸の設定 wlen=fs*0.032; %窓長 wshift=fs*0.008; % 窓シフト wnum=length(x)/wshift-3; %フレーム数の計算 wnum=floor(wnum); % 整数に直す f=figure; % 図のID取得 set(f,'doublebuffer','on');%アニメーション表示設定 M=moviein(wnum,f); %ムービーバッファ初期化 hw=hanning(wlen); for n=0:wnum-1 S=n*wshift; %フレーム開始点 xtmp=x(S+1:S+wlen); % フレームデータ取得 xtmp2=hw.*xtmp; %窓掛け Z(1:wlen)= xtmp2'; % データの転置 subplot(2,1,1);plot(t(1:S+wlen),x(1:S+wlen), ... t(S+wlen+1:end),x(S+wlen+1:end)); xlabel('time');ylabel('Amplitude'); subplot(2,1,2);plot(Z);axis([1 wlen -10^4 10^4]); xlabel('Data Points');ylabel('Amplitude'); M(:,n+1)=getframe(f);%ムービー画像の作成 disp(['n=' int2str(n)]); %フレーム番号表示 end movie(f,M);%ムービーの再生
%===================================================================
close all; % 画面クリア clear all; % データクリア [x,fs]=aiffread('kousyouAD16.aiff'); %データ読み込み t=[0:1/fs:length(x)/fs-1/fs]; %時間軸の設定 wlen=fs*0.032; %窓長 wshift=fs*0.008; % 窓シフト wnum=length(x)/wshift-3; %フレーム数の計算 wnum=floor(wnum); % 整数に直す f=figure; % 図のID取得 set(f,'doublebuffer','on');%アニメーション表示設定 M=moviein(wnum,f); %ムービーバッファ初期化 hw=hanning(wlen); fq=(1:wlen)*fs/wlen; fq=fq(1:wlen/2);%周波数軸の設定 for n=0:wnum-1 S=n*wshift; %フレーム開始点 xtmp=x(S+1:S+wlen); % フレームデータ取得 xtmp2=hw.*xtmp; y=fft(xtmp2'); Z= 20*log10(abs(y(1:wlen/2))) ; subplot(2,1,1);plot(t(1:S+wlen),x(1:S+wlen), ... t(S+wlen+1:end),x(S+wlen+1:end)); xlabel('time');ylabel('Amplitude'); subplot(2,1,2);plot(fq,Z);axis([1 fs/2 0 120]); xlabel('Frequency');ylabel('Magnitude'); M(:,n+1)=getframe(f);%ムービー画像の作成 disp(['n=' int2str(n)]); %フレーム番号表示 end movie(f,M);%ムービーの再生
%===================================================================
close all; % 画面クリア clear all; % データクリア [x,fs]=aiffread('kousyouAD16.aiff'); %データ読み込み t=[0:1/fs:length(x)/fs-1/fs]; %時間軸の設定 wlen=fs*0.032; %窓長 wshift=fs*0.008; % 窓シフト wnum=length(x)/wshift-3; %フレーム数の計算 wnum=floor(wnum); % 整数に直す f=figure; % 図のID取得 set(f,'doublebuffer','on');%アニメーション表示設定 M=moviein(wnum,f); %ムービーバッファ初期化 hw=hanning(wlen); for n=0:wnum-1 S=n*wshift; %フレーム開始点 xtmp=x(S+1:S+wlen); % フレームデータ取得 xtmp2=hw.*xtmp; y=fft(xtmp2'); Z=real(ifft(log(abs(y))));% ケプストラム subplot(2,1,1);plot(t(1:S+wlen),x(1:S+wlen), ... t(S+wlen+1:end),x(S+wlen+1:end)); xlabel('time');ylabel('Amplitude'); subplot(2,1,2);plot(Z);axis([1 wlen -1 1]); xlabel('Quefrency');ylabel('Amplitude'); M(:,n+1)=getframe(f);%ムービー画像の作成 disp(['n=' int2str(n)]); %フレーム番号表示 end movie(f,M);%ムービーの再生
%===================================================================
close all; % 画面クリア clear all; % データクリア [x,fs]=aiffread('kousyouAD16.aiff'); %データ読み込み t=[0:1/fs:length(x)/fs-1/fs]; %時間軸の設定 wlen=fs*0.032; %窓長 wshift=fs*0.008; % 窓シフト wnum=length(x)/wshift-3; %フレーム数の計算 wnum=floor(wnum); % 整数に直す f=figure; % 図のID取得 set(f,'doublebuffer','on');%アニメーション表示設定 M=moviein(wnum,f); %ムービーバッファ初期化 hw=hanning(wlen); fq=(1:wlen)*fs/wlen; fq=fq(1:wlen/2);%周波数軸の設定 for n=0:wnum-1 S=n*wshift; %フレーム開始点 xtmp=x(S+1:S+wlen); % フレームデータ取得 xtmp2=hw.*xtmp; y=fft(xtmp2'); N=30; cep=real(ifft(log(abs(y)))); lifter=[ones(1,N) zeros(1,length(cep)-2*N) ones(1,N)];%リフタの作成 lowcep=cep.*lifter;%リフタの重みつけによる低次ケプストラムの抽出 yy=fft(lowcep); Z1= 20*log10(abs(y(1:wlen/2))) ; Z2= 20*log10(exp(abs(yy(1:wlen/2))) ); subplot(2,1,1);plot(t(1:S+wlen),x(1:S+wlen), ... t(S+wlen+1:end),x(S+wlen+1:end)); xlabel('time');ylabel('Amplitude'); subplot(2,1,2);plot(fq,Z1,fq,Z2);axis([1 fs/2 0 120]); xlabel('Frequency');ylabel('Magnitude'); M(:,n+1)=getframe(f);%ムービー画像の作成 disp(['n=' int2str(n)]); %フレーム番号表示 end movie(f,M);%ムービーの再生
%===================================================================
MATLAB Getting Started with MATLAB
ユーザーズディジタル信号処理