FFT without inbuilt functions

I am trying to find fft for ECG signal. I am using .mat files. Right now I did the easy part with using inbuilt fft but I need to find fft without inbuilt funcitons beause of academic reasons. However I did not quite understand create a algorithm for fft and combine it .mat file.

clc;
load ('s0010_rem.mat')
ECGsignal = val;
Fs=1000;
L = 1000;
t = (0:length(ECGsignal)-1)/Fs;

val1 = val(1,:)/Fs; 
val2 = val(2,:)/Fs;
val3 = val(3,:)/Fs;
val4 = val(4,:)/Fs; 
val5 = val(5,:)/Fs;
val6 = val(6,:)/Fs;
val7 = val(7,:)/Fs;
val8 = val(8,:)/Fs; 
val9 = val(9,:)/Fs;
val10 = val(10,:)/Fs;

The part of fft with inbuilt:

NFFT = 2^nextpow2(L); % Next power of 2 from length of y
val1_fft = fft(val1,NFFT)/L;
val2_fft = fft(val2,NFFT)/L;
val3_fft = fft(val3,NFFT)/L;
val4_fft = fft(val4,NFFT)/L;
val5_fft = fft(val5,NFFT)/L;
val6_fft = fft(val6,NFFT)/L;
val7_fft = fft(val7,NFFT)/L;
val8_fft = fft(val8,NFFT)/L;
val9_fft = fft(val9,NFFT)/L;
val10_fft = fft(val10,NFFT)/L;
f = Fs/2*linspace(0,1,NFFT/2+1);

Answer 1

[Onur Dikilitas]

NFFT = 2^nextpow2(L); % Next power of 2 from length of y

for k=1:NFFT
    val1_fft(k) = 0;
    for n = 1:NFFT
        val1_fft(k)=val1_fft(k)+val1(n).*exp((-1j).*2.*pi.*(n-1).*(k-1)./NFFT);
    end
end
val1_fft = val1_fft/L;

[Onur Dikilitas]

I figured it out. This is a example for single one.

Answer 2

[ROBIN SINGH SOM]

% 1sep2021
% fft operation without using inbuilt function
% only capable to compute fft upto N = 128 but it can be easily expendable
clc
clear
% generating random signal of length 30
x = randi([-5,5],1,30);
% length of the signal
N = length(x);
% number of stage required
M = log2(N);
% adding zeros 
if (rem(M,1) ~= 0)
    re =rem(M,1);
    M=M-re+1;
    Ne = 2^M;
    x = [x,zeros(1,Ne-N)];
else
    Ne = N;
end
% calc. using inbuild function
x_fft =fft(x);
% bitreversing
x = bitrevorder(x);
% instialization of variables used for different stages
temp = zeros(1,Ne);
temp2 = zeros(1,Ne);
temp3 = zeros(1,Ne);
temp4 = zeros(1,Ne);
temp5 = zeros(1,Ne);
temp6 = zeros(1,Ne);
temp7 = zeros(1,Ne);
% code 
for l = 1:M
    if l==1
        for t=0:2:Ne-1
            temp(t+1:t+2) = temp(t+1:t+2) + myfun(x(t+1:t+2),2^l);
        end
        out = temp;
    end
    if l==2
        for k=0:4:Ne-1
            temp2(k+1:k+4) = temp2(k+1:k+4) + myfun(temp(k+1:k+4),2^l);
        end  
        out= temp2;
    end
    if l==3
        for k=0:8:Ne-1
            temp3(k+1:k+8) = temp3(k+1:k+8) + myfun(temp2(k+1:k+8),2^l);
        end  
        out= temp3;
    end
    if l==4
        for k=0:16:Ne-1
            temp4(k+1:k+16) = temp4(k+1:k+16) + myfun(temp3(k+1:k+16),2^l);
        end  
        out= temp4;
    end
    if l==5
        for k=0:32:Ne-1
            temp5(k+1:k+32) = temp5(k+1:k+32) + myfun(temp4(k+1:k+32),2^l);
        end  
        out= temp5;
    end
    if l==6
        for k=0:64:Ne-1
            temp6(k+1:k+64) = temp6(k+1:k+64) + myfun(temp5(k+1:k+64),2^l);
        end  
        out= temp6;
    end
    if l==7
        for k=0:128:Ne-1
            temp7(k+1:k+128) = temp7(k+1:k+128) + myfun(temp6(k+1:k+128),2^l);
        end  
        out= temp7;
    end
end
x % values of x 

x = 1×32

-4 0 2 2 0 5 -4 -5 4 1 1 -3 -1 1 -5 0 1 -4 2 -5 -1 4 4 -3 -2 1 -3 5 -5 5

x_fft % fft calc. using inbuild function

x_fft = 
  -5.0000 + 0.0000i -10.1699 +11.9096i  -6.6514 -20.7960i  19.7950 - 2.3288i  -2.3640 - 0.2218i -10.4373 -17.8994i  -4.0294 + 6.9798i  11.9712 - 0.8133i  -2.0000 + 5.0000i -38.6463 + 7.0083i -14.7990 - 4.0496i  -6.2870 -17.4688i  10.3640 +15.7782i  -2.1003 -13.2119i  -6.5201 +24.1745i   3.8745 + 8.4175i  -7.0000 + 0.0000i   3.8745 - 8.4175i  -6.5201 -24.1745i  -2.1003 +13.2119i  10.3640 -15.7782i  -6.2870 +17.4688i -14.7990 + 4.0496i -38.6463 - 7.0083i  -2.0000 - 5.0000i  11.9712 + 0.8133i  -4.0294 - 6.9798i -10.4373 +17.8994i  -2.3640 + 0.2218i  19.7950 + 2.3288i  -6.6514 +20.7960i -10.1699 -11.9096i

out   % fft calc. using myfunction (DIT)

out = 
  -5.0000 + 0.0000i -10.1699 +11.9096i  -6.6514 -20.7960i  19.7950 - 2.3288i  -2.3640 - 0.2218i -10.4373 -17.8994i  -4.0294 + 6.9798i  11.9712 - 0.8133i  -2.0000 + 5.0000i -38.6463 + 7.0083i -14.7990 - 4.0496i  -6.2870 -17.4688i  10.3640 +15.7782i  -2.1003 -13.2119i  -6.5201 +24.1745i   3.8745 + 8.4175i  -7.0000 - 0.0000i   3.8745 - 8.4175i  -6.5201 -24.1745i  -2.1003 +13.2119i  10.3640 -15.7782i  -6.2870 +17.4688i -14.7990 + 4.0496i -38.6463 - 7.0083i  -2.0000 - 5.0000i  11.9712 + 0.8133i  -4.0294 - 6.9798i -10.4373 +17.8994i  -2.3640 + 0.2218i  19.7950 + 2.3288i  -6.6514 +20.7960i -10.1699 -11.9096i

figure()
hold on
stem(abs(out),'filled','LineStyle',"--","Marker","diamond","Color",'b',"LineWidth",1)
title("Without Inbuilt Function","FontSize",15)
hold off

figure()
hold on
stem(abs(x_fft),'filled','LineStyle',"-.",'Marker',"*","color",'r',"LineWidth",1)
title("With Inbuilt Function","FontSize",15)
hold off

figure()
hold on
stem(abs(x_fft),'filled','LineStyle',"-.",'Marker',"*","color",'r',"LineWidth",1)
stem(abs(out),'filled','LineStyle',"none","Marker","diamond","Color",'b',"LineWidth",1)
title("Overlapping","FontSize",15)
hold off

   

function out = myfun(inp,N)
    twi = twiddle(N);
    inp1 = [inp(1:N/2),inp(1:N/2)];
    inp2 = [inp(N/2+1:N),inp(N/2+1:N)];
    out = zeros(1,N);
    for i=1:N
        out(i) = inp1(i) + twi(i)*inp2(i);
    end
end

function out = twiddle(N)
out = zeros(1,N);
for k=1:N
  out(k) = exp(-1j*2*pi*(k-1)/N);
end
end



% author: Robin Singh Som 

Answer 3

[ROBIN SINGH SOM]

% 7sep2021
% fft operation without using inbuilt function
clc
clear
q = 10000;
% generating sin signale with 10000 samples
x = sin(2*pi*0.2*(0:q-1));
% length of the signal
N = length(x);
M = log2(N);
% adding zeros 
if (rem(M,1) ~=0)
    re = rem(M,1);
    M = M-re+1;
    Ne = 2^M;
    x = [x,zeros(1,Ne-N)];
else
    Ne = N;
end
N = Ne;
% number of stages
M = log2(N);
% calc. using inbuild function
x_fft =fft(x);
% bitreversing
x = bitrevorder(x);
% code 

for l =1:M
    k = 2^l;
    w = exp((-1j*2*pi*(0:k-1))./(k));
         for t=0:k:N-1 
             z(t+1:t+k) = [x(t+1:t+k/2)+x(t+k/2+1:t+k) .* w(1:k/2) , x(t+1:t+k/2)+x(t+k/2+1:t+k).*w(k/2+1:k)] ;            
         end
         x = z;
         z = 0;         
 end
out = x;
figure()
hold on
stem((0:N-1/N),(abs(out)),'filled','LineStyle',"--","Marker","diamond","Color",'b',"LineWidth",1)
title("Without Inbuilt Function","FontSize",15)
hold off

figure()
hold on
stem((0:N-1/N),(abs(x_fft)),'filled','LineStyle',"-.",'Marker',"*","color",'r',"LineWidth",1)
title("With Inbuilt Function","FontSize",15)
hold off

figure()
hold on
stem(abs(x_fft),'filled','LineStyle',"-.",'Marker',"*","color",'r',"LineWidth",1)
stem(abs(out),'filled','LineStyle',"none","Marker","diamond","Color",'b',"LineWidth",1)
title("Overlapping","FontSize",15)
hold off

%robin singh