Need FFT Code for Matlab (not built in)

Question

John 2013-3-15

1
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/67283-need-fft-code-for-matlab-not-built-in

评论： Marco Marcon 2024-4-28

Does anyone have FFT code for Matlab? I don't want to use the built-in Matlab function.

5 个评论
显示 3更早的评论隐藏 3更早的评论

John 2013-3-15

Okay I am looking for the O(n log n) algorithm for DFT then =)

Yasin Arafat Shuva 2021-7-19

编辑：Rik 2021-7-19

在 MATLAB Online 中打开

clear variables
close all
%%Implementing the DFT in matrix notation
xn = [0 2 4 6 8];        % sequence of x(n)
N = length(xn);        % The fundamental period of the DFT
n = 0:1:N-1;           % row vector for n
k = 0:1:N-1;           % row vecor for k
WN = exp(-1i*2*pi/N);  % Wn factor 
nk = n'*k;             % creates a N by N matrix of nk values
WNnk = WN .^ nk;       % DFT matrix
Xk = xn * WNnk;        % row vector for DFT coefficients
disp('The DFT of x(n) is Xk = ');
disp(Xk)
magXk = abs(Xk);       % The magnitude of the DFT
%%Implementing the inverse DFT (IDFT) in matrix notation
n = 0:1:N-1;           
k = 0:1:N-1;           
WN = exp(-1i*2*pi/N);  
nk = n'*k;             
WNnk = WN .^ (-nk);    % IDFS matrix
x_hat = (Xk * WNnk)/N; % row vector for IDFS values
disp('and the IDFT of x(n) is x_hat(n) = ');
disp(real(x_hat))
% The input sequence, x(n)
subplot(3,1,1); 
stem(k,xn,'linewidth',2); 
xlabel('n'); 
ylabel('x(n)')
title('plot of the sequence x(n)')
grid
% The DFT
subplot(3,1,2); 
stem(k,magXk,'linewidth',2,'color','r'); 
xlabel('k'); 
ylabel('Xk(k)')
title('DFT, Xk')
grid
% The IDFT
subplot(3,1,3); 
stem(k,real(x_hat),'linewidth',2,'color','m'); 
xlabel('n'); 
ylabel('x(n)')
title('Plot of the inverse DFT, x_{hat}(n) = x(n)')
grid

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Youssef Khmou 2013-3-15

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/67283-need-fft-code-for-matlab-not-built-in#answer_78731

在 MATLAB Online 中打开

hi John

Yes there are many versions of Discrete Fourier Transform :

 % function z=Fast_Fourier_Transform(x,nfft)
%
% N=length(x);
% z=zeros(1,nfft);
% Sum=0;
% for k=1:nfft
%     for jj=1:N
%         Sum=Sum+x(jj)*exp(-2*pi*j*(jj-1)*(k-1)/nfft);
%     end
% z(k)=Sum;
% Sum=0;% Reset
% end
% return

this is a part of the code , try my submission it contains more details :

http://www.mathworks.com/matlabcentral/fileexchange/40002-psd-power-spectral-density-and-amplitude-spectrum-with-adjusted-fft

4 个评论
显示 2更早的评论隐藏 2更早的评论

Abdullah Khan 2013-11-28

编辑：Abdullah Khan 2013-11-28

在 MATLAB Online 中打开

This is also DFT

The complexity of the code is still N x N

because

% nfft=2^ceil(log2(N)) = nfft (on upper DTFT code )

or

% 2^ceil(log2(N))= N

Walter Roberson 2017-9-6

"FFT" stands for "Fast Fourier Transform", which is a particular algorithm to make Discrete Fourier Transform (DFT) more efficient.

请先登录，再进行评论。

Answer 2

Tobias 2014-2-11

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/67283-need-fft-code-for-matlab-not-built-in#answer_123974

编辑：Walter Roberson 2017-9-6

在 MATLAB Online 中打开

For anyone searching an educating matlab implementation, here is what I scribbled together just now:

function X = myFFT(x)
%only works if N = 2^k
N = numel(x);
xp = x(1:2:end);
xpp = x(2:2:end);
if N>=8
    Xp = myFFT(xp);
    Xpp = myFFT(xpp);
    X = zeros(N,1);
    Wn = exp(-1i*2*pi*((0:N/2-1)')/N);
    tmp = Wn .* Xpp;
    X = [(Xp + tmp);(Xp -tmp)];
else
    switch N
        case 2
            X = [1 1;1 -1]*x;
        case 4
            X = [1 0 1 0; 0 1 0 -1i; 1 0 -1 0;0 1 0 1i]*[1 0 1 0;1 0 -1 0;0 1 0 1;0 1 0 -1]*x;
        otherwise
            error('N not correct.');
    end
end
end

4 个评论
显示 2更早的评论隐藏 2更早的评论

Siddhant Sharma 2020-11-2

can anyone explain X equal to (Xp-tmp)?

Marco Marcon 2024-4-28

在 MATLAB Online 中打开

The total length of the final X will be N.

The first

values will be given by

+

(remember that

,

and

have length

elements and both

,

are periodic of

); it can be written as

for

, where

.

While the last

values will be given by

for

.

Since

we can get the second part of X(last

values) as:

that is the reason for adding tmp in the first vector part and subtracting it in the second part.

For more details: Wikipedia: Cooley-Tukey algorithm

请先登录，再进行评论。

Answer 3

Jan Afridi 2017-9-6

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/67283-need-fft-code-for-matlab-not-built-in#answer_280448

编辑：Jan Afridi 2017-11-6

I think I can answer your question. Check the given link may be it will help you. And if you found any error then please inform me... check the link Matlab code for Radix 2 fft

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Jan Suchánek 2018-12-6

Hello i would like to ask you what had to be different if you would like to use it for harmonical function.

I tried to implement your solution but it was working only for random vectors but not for any type of harmonical.

i was using this for the harmonical function.

______________________

t2=0:1:255;

h1=sin(t2);

x=h1;

X11 = myFFT(x);

______________________

it gave me this error:

Error using *

Inner matrix dimensions must agree.

Error in myFFT (line 18)

X = [1 0 1 0; 0 1 0 -1i; 1 0 -1 0;0 1 0 1i]*[1 0 1 0;1

0 -1 0;0 1 0 1;0 1 0 -1]*x;

Error in myFFT (line 7)

Xp = myFFT(xp);

Error in myFFT (line 7)

Xp = myFFT(xp);

Error in myFFT (line 7)

Xp = myFFT(xp);

Error in myFFT (line 7)

Xp = myFFT(xp);

Error in myFFT (line 7)

Xp = myFFT(xp);

Error in myFFT (line 7)

Xp = myFFT(xp);

Error in ZBS_projekt (line 246)

X11 = myFFT(x);

请先登录，再进行评论。

Answer 4

Ryan Black 2019-3-11

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/67283-need-fft-code-for-matlab-not-built-in#answer_364830

编辑：Ryan Black 2020-8-26

Function I wrote that optimizes the DFT or iDFT for BOTH prime and composite signals...

https://www.mathworks.com/matlabcentral/fileexchange/70563-rader-s-mixed-radix-fft

The forward transform is triggered by -1 and takes the time-signal as a row vector:

[Y] = fftmodule(y,-1)

The inverse transform auto-normalizes by N, is triggered by 1 and takes the frequency-signal as a row vector:

[y] = fftmodule(Y,1)

Methodology:

The algorithm decimates to N's prime factorization following the branches and nodes of a factor tree. In formal literature this may be referred to as Mixed Radix FFT, but its really just recursive decimation of additive groups and this method is easily derivable via circular convolutions :)

At the prime tree level, algorithm either performs a naive DFT or if needed performs a single Rader's Algorithm Decomposition to (M-1), zero-pads to power-of-2, then proceeds to Rader's Convolution routine (not easy to derive). Finally it upsamples through the origianlly strucutured nodes and branches incorporating twiddle factors for the solution.