初めまして
私は今、matlabで離散フーリエ変換をfft関数などを使わずにプログラムを書いています。「dft_size×dft_size の行列」と「dft_size×1 の行列」の積により DFT 計算を行うプログラムを書きたいのですが上手く書けなくて困っています。

下記プログラムの
X_real=zeros(dft_size,1);
X_imag=zeros(dft_size,1);のX_real=zeros()内を(dft_size,dft_size)に変えて、
scriptファイルの
A=zeros(dft_size,1);
T=zeros(dft_size,1);のA()内を(dft_size,dft_size)に書き換えればよいと思ったのですが、上手くいきませんでした。

アドバイスを頂けたら幸いです。

以下、プログラムリスト
//////

DFT.m

function [X_real,X_imag]=DFT(x_real,x_imag,dft_size)
X_real=zeros(dft_size,1);
X_imag=zeros(dft_size,1);
for k=1:dft_size,
for n=1:dft_size,
w_real=cos(-2pi(k-1)(n-1)/dft_size);
w_imag=sin(-2pi*(k-1)(n-1)/dft_size);
X_real(k,1)=X_real(k,1)+w_realx_real(n,1)-w_imagx_imag(n,1);
X_imag(k,1)=X_imag(k,1)+w_realx_imag(n,1)+w_imag*x_real(n,1);
end
end
////////

script.m

clear;
dft_size=60;
x_real=audioread('test.wav',[1 dft_size]);
x_imag=zeros(dft_size,1);
[X_real,X_imag]=DFT(x_real,x_imag,dft_size);
A=zeros(dft_size,1);
T=zeros(dft_size,1);
frequency=zeros(dft_size,1);
for k=1:dft_size,
A(k,1)=sqrt(X_real(k,1)*X_real(k,1)+X_imag(k,1)*X_imag(k,1));
T(k,1)=atan(X_imag(k,1)/X_real(k,1));
frequency(k,1)=k-1;
end
subplot(2,1,1);
stem(frequency,A);
subplot(2,1,2);
stem(frequency,T);