FFTのコマンドを用いない、DFTのアルゴリズムでのプログラミングについて
32 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
初めまして。
私は、現在MATLAB上でFFTコマンドを用いずにDFTのプログラムを作成中です。
ある20000個のデータを解析したいのですが、プログラムがうまく実行できません。
以下が作成中のプログラムとなります。
使用しているDFTの定義式は、
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/205705/image.png)
です。
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のプログラム技術や知識はあまりないので、丁寧に教えてくださると助かります。
1 commentaire
Réponse acceptée
maeda
le 26 Fév 2019
質問にあるDFTの定義式に従ってサンプルプログラムを作りました。お試しください。
% 最初に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)')
1 commentaire
maeda
le 26 Fév 2019
Modifié(e) : maeda
le 7 Mar 2019
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)')
Plus de réponses (0)
Voir également
Catégories
En savoir plus sur フーリエ解析とフィルター処理 dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!