FFTのコマンドを用​いない、DFTのアル​ゴリズムでのプログラ​ミングについて

初めまして。

私は、現在MATLAB上でFFTコマンドを用いずにDFTのプログラムを作成中です。

ある20000個のデータを解析したいのですが、プログラムがうまく実行できません。

以下が作成中のプログラムとなります。

使用しているDFTの定義式は、

　　　（k=0,1,2,…,N-1）

です。

clear;
close all;

filename1 ='C.csv';
sheet=1;
xlRange='A:A';
A=xlsread(filename1,sheet,xlRange);

N=20000;

for k=0:1:N-1
    syms I
    S1=symsum(A(I+1,1)*exp(-1j*(2*pi/N)*k*I),I,0,N-1);
    ck=(1/N)*S1;
    rck=real(ck);
end

aCk=abs(rck);
plot(aCk)

このプログラムを実行した結果、

エラー: sym/subsindex (line 769)
インデックス付けまたは関数定義が無効です。関数を定義する場合は、必ず引数をシンボリッ
ク変数にして、関数の本体を SYM 式にしてください。インデックス付けの場合、入力は数
値、論理値または ':' でなければなりません。

エラー: DFT (line 16)
    S1=symsum(A(I+1,1)*exp(-1j*(2*pi/N)*k*I),I,0,N-1);

と表示されました。

DFTを行い、振幅スペクトルおよび位相スペクトルをグラフで表示するにはどうしたらよいでしょうか。

皆様のお力を貸していただければと思います。

なお、MATLABのプログラム技術や知識はあまりないので、丁寧に教えてくださると助かります。

[michio]

質問の投稿、ありがとうございました。

回答の内容で課題や疑問が解決されましたら、

ぜひ「この回答を採用」ボタンのクリックをお願いいたします。

Answer 1

[maeda]

質問にあるＤＦＴの定義式に従ってサンプルプログラムを作りました。お試しください。

% 最初にDFTを実行したい信号yを準備します
% 0.2 Hz、振幅 1 の正弦波に30%のランダムノイズを載せた16秒間の波を準備します
Fs = 1; % サンプリング周波数は1Hz
t = (0:Fs:15)'; % 16秒間の時間ベクトルを準備
y = sin(2*pi*0.2*t) + rand(16,1)*0.3; % 0.2Hzの波を準備
N = length(y); % 信号の長さ
% DFTをかけたい信号を表示します
figure;plot(t,y);grid

% DFTを定義式に従って計算します
i = (0 : N-1);         % yのデータの個数に対応するiを準備
k = (0 : N-1)';        % 周波数の個数に対応するkを準備
expart = exp(complex(0,-(2*pi/N)*k*i)); %　exp(-j*(2*pi/N)*k*i)に対応
C = 1/N * sum(y.*expart);             % 離散フーリエ変換を実行します

% DFT後の結果をプロットします
f = Fs*(0:N-1)/N;
figure;plot(f,C);grid
xlabel('frequency(Hz)')
ylabel('abs(C)')

[maeda]

16個のデータを処理する為にすべて行列で解いたのが先ほどのスクリプトです。20000個のデータを処理するためには、20000×20000の行列を一つ作ることになり大容量のメモリが必要になります。for文を一回使ってメモリの節約を試みました。

Fs = 1; % サンプリング周波数は1Hz
L = 20000; % データの個数
t = (0:Fs:L-1)'; % 16秒間の時間ベクトルを準備
y = sin(2*pi*0.2*t) + rand(L,1)*0.3; % 0.2Hzの波を準備

% DFTをかけたい信号を表示します
figure;plot(t,y);grid

% DFTがかけられるようにyのデータ個数を2のべき乗になるようにする。yの末尾にゼロを加える
n = 2^nextpow2(L);    %DFTがかけられるようにするためのデータの個数2000なら2048個になります。
y = vertcat(y,zeros(n-L,1));
N = length(y);

% DFTを式に従って計算します
i = (0 : N-1);         % yのデータの個数に対応するiを準備
C = zeros(N,1);
for k = 0 : 1: N-1
  expart = exp(complex(0,-(2*pi/N)*k*i))'; %　exp(-j*(2*pi/N)*k*i)に対応
  C(k+1,1) = 1/N * sum(y.*expart);             % 離散フーリエ変換を実行します
end
% DFT後の結果をプロットします
f = Fs*(0:N-1)/N;
figure;plot(f,abs(C));grid
xlabel('frequency(Hz)')
ylabel('abs(C)')