課題が終わっていない件

kazuma.です。

今日は30分チャレンジやりたかったのですが、できてません。
課題を先に終わらせておきます。（締切今日の12時ヤバい・・・。）

朝が弱くて今日の朝8時に起きたのに何もできず・・・。
朝大切にしないとですね。

ちなみに今日はscilabでコードを書いてました。なんかC言語に似てます。
でも人工知能的なことあんまり好きじゃないです・・・。（笑）

clear
FS=10000;
fft_len=512;               // FFTのサンプル数
start=5000;                // 分析位置
len=300;                   // 分析長
order=20;                  // 分析次数
pre_emp=0.0;               // プリエンファシス
print(%io(2),'start= '); start   = read(%io(1),1,1);
print(%io(2),start,len,fft_len,order); 

  x=loadwave('3.wav');   // 1.wavというファイルの読み込み

// Hamming窓掛け&プリエンファシス
for i=1:len 
    win(i) = 0.54 - 0.46 * cos(2 * %pi * i / len);
     x1(i) = (x(i+start)-pre_emp*x(i-1+start)) * win(i);
end
for i=len+1:fft_len   // ゼロ詰め 
     x1(i) = 0;
    win(i) = 0;
end

// DFT対数スペクトルの算出
fft_spc=20*log10(abs(fft(x1,-1)));

// 自己相関関数の算出
j=1;
for i=0:order
    r(j)=0;    
  for n=1:len-i
    r(j)=r(j)+x1(n)*x1(n+i);
  end
   j=j+1;
end

// Levinson算法によるLPC分析
[ar,sigma2,rc]=lev(r); // sigma2=r(0)+_sum_{i=1}^order(r(i)*a(i))

// LPC対数スペクトルの算出
       a(1)=1;
     for i=1:order
       a(i+1)=ar(i);
     end
     for i=order+1:fft_len-1  // ゼロ詰め 
       a(i+1)=0;
     end
    ar_spc=-20*log10(abs(fft(a,-1)))+10*log10(sigma2);

// 残差のスペクトルの算出
   for n=1:len
       res(n)=x(start+n);
     for i=1:order
       res(n)=res(n)+ar(i)*x(start+n-i);
     end
   end
   for n=len+1:fft_len
       res(n)=0;
   end
       res_spc=20*log10(abs(fft(res,-1)));

// スペクトルの描画
 xset('window',1);  clf(1); 
 tics=[2,4,2,4];
 Hz_flag=1;
if Hz_flag == 0
 rect=[1,min(fft_spc),fft_len/2,max(fft_spc)];
 plotframe(rect,tics,[%f,%f],['LPC','Freq.','Amp.[dB]'],[0,0,1.0,0.5]);
 n=1:fft_len/2;
 plot2d(n,fft_spc(n),1,"000");
 plot2d(n, ar_spc(n),2,"000");
 plot2d(n,res_spc(n),3,"000");
else
 rect=[0,min(min(fft_spc),min(ar_spc),min(res_spc)),FS/2,max(max(fft_spc),max(ar_spc),max(res_spc))];
 plotframe(rect,tics,[%f,%f],['LPC','Freq.[Hz]','Amp.[dB]'],[0,0,1.0,0.5]);
 Hz=(0:fft_len/2)/(fft_len/2)*(FS/2); // 修正
 fft_spc=fft_spc(1:(fft_len/2+1));
 ar_spc = ar_spc(1:(fft_len/2+1));
 res_spc=res_spc(1:(fft_len/2+1));
 plot2d(Hz,fft_spc,1,"000");
 plot2d(Hz, ar_spc,2,"000");
 plot2d(Hz,res_spc,3,"000");
// set(gca(),'data_bounds', [0,min(min(fft_spc),min(ar_spc),min(res_spc));FS/2,max(max(fft_spc),max(ar_spc),max(res_spc))]); // 追加
// xgrid();
end
// 極を算出し描画
 HAR=poly(a(1:order+1),'z','coeff');
 pp=roots(HAR);
 for i=1:order
  pp(i)=1/pp(i);
 end
 x=0:0.1:2*%pi;
// xbasc();
 rect=[-1,-1,1,1];
 tics=[2,5,2,5];
 plotframe(rect,tics,[%f,%f],["Unit Circle","Re.","Im."],[0.25,0.5,0.5,0.5]);
 plot2d(cos(x),sin(x),1,"000");
 ra=real(pp); ia=imag(pp);
 plot2d(ra,ia,-3);
 xgrid();

||