高速フーリエ変換

信号

FFT

clear; close all; clc;

N=2^8;
fs = 100;
T=1/fs;
f1=10;f2=20;f3=40;
A1=1;A2=0.5;A3=0.8;
t=0:T:(N-1)*T;
scale = 0.3;
noise = rand(1,length(t))*scale;
x=A1*sin(2*pi*f1*t)+A2*sin(2*pi*f2*t)+A3*sin(2*pi*f3*t)+noise;
plot(t,x);
xlim([0,0.5]);

F = myFFT(x);
amp = 2*abs(F)/N;
freq = linspace(0,fs,N);
figure();
plot(freq(1:N/2+1),amp(1:N/2+1));

function hatx=myFFT(x)
  N = length(x);
  hatx = zeros(1,N);
  k = 1:floor(N/2);
  W = exp(-1j*2*pi*k/N);
  if N==2
      hatx(1) = x(1) + x(2);
      hatx(2) = x(1) - x(2);
  end
  if N >=4      
     L = myFFT(x(1:2:end));
     R = myFFT(x(2:2:end));
     hatx(1:floor(N/2))   = L + W.*R;
     hatx(floor(N/2)+1:N) = L - W.*R;
  end   
end