ですから、ｋ＝8の場合、N=8です。
これは、行列で表現すると、以下の計算になります。
ｘ（ｋ）は、入力。
X(k)は、出力です。(「ディジタルサウンド処理入門青木先生CQ出版」）
=
注意する点は、
X(k)のKは、直接に、周波数を示すものでは、無くて
いま、sampling周期を、Tsとすると、
N(今は、N=k=8）を使って
１番長い周期N*Ts（これは、周波数で言えば、１番低い周波数FF=1/(N*Ts) ）
から、整数倍の、k番目の周波数
F*ｋ[ｋ＝０の時、直流]
を、表すことです。
入力から、ｋ個の点を取り出す、ちゅうことは、
そこで、信号を打ち切って、その周期を、最長の周期として、その整数倍の周波数を計算する
ことなんやろかな....
sampling周期Tsのk倍（ｋ*Ts）を、最長の周期（つまり、基本の周波数F=1/(k*Ts））
として、その整数倍の周波数F*kの成分を調べる、ことに、なります。
ですので
k=1：基本周波数の成分k=2：基本周波数の2倍の成分k=3：基本周波数の3倍の成分.......
各成分を計算することに、なります。(k=0は、直流成分のようですね）
「標本化周波数ｆｃをNで等分割した周波数に限り、値をもつことが知られています。」（ディジタルサウンド入門青木先生ｐ.30）
これは、知らんかったな...
すると、基本周波数の成分X(1)は
で、表されることに、なります。
と、書いたものの、実は、何故そうなのか？まだ、解らない.....
井澤先生のHPを拝見しながら、考えていますが、まだ、、意味合いが、解らん.....涙 (-_-)
ほな、基本周波数の2倍の成分X(2)は
や、さかい
あかん、....解らん.....
複素数のままで、計算します。
まず、回転の行列を作る関数
function w=make_W(n)
w=zeros(n,n);
for k=1:n
for l=1:n
w(k,l)=exp(-i*2*pi/n)^((k-1)*(l-1));
end
end
endfunction
次に、信号も複素数のままで、行列の計算をして、DFTする。
function dft_out=mydft_(x)
n=length(x);
w=zeros(n,n);
dft_out=zeros(n,1);
w=make_W(n);%さっき作った関数を呼び出す
dft_out=w* (x)';
end
endfunction
そんで、これをDFTして、結果を、グラフ表示するmファイル
clear
x_real=wavread('test.wav');
x=x_real(1:60); %1から60までの点をxに代入
z=mydft_(x');
for k=1:60
z_abs(k)=abs(z(k));%振幅を計算
z_atan(k)=atan(imag(z(k))/real(z(k))); %位相を計算
freq_(k)=k-1;%freq_は１から始まるのを、0からに修正
end
subplot(2,1,1);
stem(freq_,z_abs');
subplot(2,1,2);
stem(freq_,z_atan');
このmydft_が、うまく動作する事を、既存のDFT_関数で、確かめておく。
下は、make_w関数を確かめています。
只、インストールした、だけでは、文字が小さすぎます。
これを、修正して、文字を大きくするには
Octaveのショートカットのプロパティを出して（右クリックして）
リンク先と書いた所の、
Console-12
を
Console-16
に変えておきます。（16くらいで、丁度ええんちゃう？）
ねっ、見やすくなったでしょう。
そうそう、gnuplotの画面の、左上、アイコンをクリックすると、
グラフの色、フォント、線幅等を、変えて、defaultとして、保存できます。
Scilab、Octaveも、ええけど、 Eulerも、素晴らしいですね....
Eulerが、私の好みです。（また、へそ曲がり、ばっかり言うて...）
みなさん、ありがとうございます。
H.18.4.18