Contents
FFT algoritm
clear; % clears all previus values from memory +clc; % clear command window +fs = 44100; % samplinf freq. +fftLength=32; % windowlength + +% signal frequencies +data_length = 8; % data length in FPGA calculations +max = 2^(data_length-1) - 1 ; % max aplitude 2^n /2 + +f1 = 1000; +a1 = max/2; + +f2 = 0; +a2 = max/4; + +f3 = 8000; +a3 = max/2; + +% calculating signals +comp1 = a1 * cos(2*pi*f1*[0:1/fs:1]); +comp2 = a2 * cos(2*pi*f2*[0:1/fs:1]); +comp3 = a3 * cos(2*pi*f3*[0:1/fs:1]); + +% calculatin vector values for step function +d1 = ones(1, 24); +d2 = 0.*ones(1, 1000 ); + +%data = [ d1 , d2]; % creates vector with step function +data = comp1 + comp2 + comp3; % creates vector from 3 sin functions + +figure(1) % plots separete sin functions +plot ( comp1, '-'); +hold on; +plot ( comp2, '-'); +plot ( comp3, '-'); +xlim([1 50]) +title('Separete SIN functions') +ylabel('magnitude'), xlabel('time') +hold off; + +figure(2) % plots signal for fft +plot ( data); +title('Signal for FFT analysis FFT') +ylabel('magnitude'), xlabel('time') +xlim([1 100]) + +figure(3) % plots resultinf fft from Matlab functions +ft =fft(data,fftLength); +ftMag=abs(ft(1:fftLength/2)); +stem (ftMag) +title('Linear Magnitude FFT') +ylabel('magnitude'), xlabel('kHz') + +xt = xticks; % returns the current x-axis tick values as a vector +fstep = fs/fftLength; % tick of f axis in f domain +xtnew = round((xt-1)*fstep/1000, 1) ; % calculate new tick in kHz +xticklabels(xtnew) % set new tick labels + +% figure(4) % plots resultinf fft(in dB) from Matlab functions +% ft =fft(data,fftLength); +% ftMag=abs(ft(1:fftLength/2)); +% plot (20*log10(ftMag)) +% title('dB Magnitude') +% ylabel('dB'), xlabel('kHz') +% +% xt = xticks; % returns the current x-axis tick values as a vector +% fstep = fs/fftLength; % tick of f axis in f domain +% xtnew = round((xt-1)*fstep/1000, 1) ; % calculate new tick in kHz +% xticklabels(xtnew) % set new tick labels +
Data preparation for FFT
% reverse bit calulation +bits = length(dec2bin( fftLength - 1 )); % how many bits in binary number +rev_bit_dec = zeros(1,fftLength); % create vektor size of fftlength + +for n=1:fftLength + bin_num = dec2bin(n-1 , bits); % converting to binary number + rev_bit = []; % create empty vector + for k=bits:-1:1 + rev_bit = [rev_bit , bin_num(k)]; + end + rev_bit_dec(n) = bin2dec(rev_bit) ; % add 1 to match Matlab numbering +end + +% creating array +% create empty array to store values in reverse bit order +stage = zeros(bits + 1,fftLength); + +for n=1:fftLength + stage(1,n) = data(rev_bit_dec(n)+1); +end + +% Calculating W twiddling factor for all stages +for n = 1 : fftLength/2 + W(n) = exp(-1i * (n-1) * 2 * pi/ fftLength ); +end + +% convert to fixed point mumber -> sfi(v,w,f) returns a signed fixed-point object with value v, word length w, and fraction length f. +Wr = sfi(real(W),data_length,data_length-2); +Wi = sfi(imag(W),data_length,data_length-2); + +st_real = sfi(real(stage) , data_length + 3 , 0); +st_imag = sfi(imag(stage) , data_length + 3 , 0); +% temp values for multiplaying with W twiddling factor +st_real_tmp = sfi(real(zeros(bits + 1,fftLength)) , data_length + 3 , 0); +st_imag_tmp = sfi(imag(zeros(bits + 1,fftLength)) , data_length + 3 , 0); +
First stage
for n = 1 : 2^1 : fftLength + % Even + stage(2,n) = stage(1,n) + stage(1,n+1); + % Odd + stage(2,n+1) = stage(1,n) - stage(1,n+1); +end + +% calculations using separate real and imaginary numbers +for n = 1 : 2^1 : fftLength + % Even + st_real(2,n) = st_real(1,n) + st_real(1,n+1); + % imag is 0 + % Odd + st_real(2,n+1) = st_real(1,n) - st_real(1,n+1); + % imag is 0 +end +
Second stage
% Calculating W twiddling factor +for n = 1 : 2 + Wn(n) = exp(-1i * (n-1) * 2 * pi/ 4 ); +end + +% calculate next stage values +for n = 1 : 2^2 : fftLength + % Even pair + stage(3,n+0) = stage(2,n+0) + Wn(1)*stage(2,n+2); + stage(3,n+1) = stage(2,n+1) + Wn(2)*stage(2,n+3); + % Odd par + stage(3,n+2) = stage(2,n+0) - Wn(1)*stage(2,n+2); + stage(3,n+3) = stage(2,n+1) - Wn(2)*stage(2,n+3); +end + +% calculations using separate real and imaginary numbers +for n = 1 : 2^2 : fftLength + % Even pair + st_real(3,n+0) = st_real(2,n+0) + st_real(2,n+2); + % imag is 0 + st_real(3,n+1) = st_real(2,n+1) ; % real is 0 + st_imag(3,n+1) = -1 * st_real(2,n+3); % mult -j + % Odd par + st_real(3,n+2) = st_real(2,n+0) - st_real(2,n+2); + % imag is 0 + st_real(3,n+3) = st_real(2,n+1) ; % real is 0 + st_imag(3,n+3) = st_real(2,n+3); % mult -j +end +
Therd stage
% Calculating W twiddling factor +for n = 1 : 4 + Wn(n) = exp(-1i * (n-1) * 2 * pi/ 8 ); +end + +% calculate next stage values +for n = 1 : 2^3 : fftLength + for k = 0 : 3 + % Even pair + stage(4,n+k) = stage(3,n+k) + Wn(k+1)*stage(3,n+k+4); + % Odd par + stage(4,n+k+4) = stage(3,n+k) - Wn(k+1)*stage(3,n+k+4); + end +end + +% calculations using separate real and imaginary numbers +for n = 1 : 2^3 : fftLength + for k = 0 : 3 + st_real_tmp(3,n+k+4) = ( Wr(k*4+1) * st_real(3,n+k+4) ) - ( Wi(k*4+1) * st_imag(3,n+k+4) ); + st_imag_tmp(3,n+k+4) = ( Wi(k*4+1) * st_real(3,n+k+4) ) + ( Wr(k*4+1) * st_imag(3,n+k+4) ); + end +end +for n = 1 : 2^3 : fftLength + for k = 0 : 3 + % Even pair + st_real(4,n+k) = st_real(3,n+k) + st_real_tmp(3,n+k+4); + st_imag(4,n+k) = st_imag(3,n+k) + st_imag_tmp(3,n+k+4); + % Odd par + st_real(4,n+k+4) = st_real(3,n+k) - st_real_tmp(3,n+k+4); + st_imag(4,n+k+4) = st_imag(3,n+k) - st_imag_tmp(3,n+k+4); + end +end +
4th stage
% Calculating W twiddling factor +for n = 1 : 8 + Wn(n) = exp(-1i * (n-1) * 2 * pi/ 16 ); +end + +% calculate next stage values +for n = 1 : 2^4 : fftLength + for k = 0 : 7 + % Even pair + stage(5,n+k) = stage(4,n+k) + Wn(k+1)*stage(4,n+k+8); + % Odd par + stage(5,n+k+8) = stage(4,n+k) - Wn(k+1)*stage(4,n+k+8); + end +end + +% calculations using separate real and imaginary numbers +for n = 1 : 2^4 : fftLength + for k = 0 : 7 + st_real_tmp(4,n+k+8) = ( Wr(k*2+1) * st_real(4,n+k+8) ) - ( Wi(k*2+1) * st_imag(4,n+k+8) ); + st_imag_tmp(4,n+k+8) = ( Wi(k*2+1) * st_real(4,n+k+8) ) + ( Wr(k*2+1) * st_imag(4,n+k+8) ); + end +end +for n = 1 : 2^4 : fftLength + for k = 0 : 7 + % Even pair + st_real(5,n+k) = st_real(4,n+k) + st_real_tmp(4,n+k+8); + st_imag(5,n+k) = st_imag(4,n+k) + st_imag_tmp(4,n+k+8); + % Odd par + st_real(5,n+k+8) = st_real(4,n+k) - st_real_tmp(4,n+k+8); + st_imag(5,n+k+8) = st_imag(4,n+k) - st_imag_tmp(4,n+k+8); + end +end +
5th stage
% Calculating W twiddling factor +for n = 1 : 16 + Wn(n) = exp(-1i * (n-1) * 2 * pi/ 32 ); +end + +% calculate next stage values +for n = 1 : 2^5 : fftLength + for k = 0 : 15 + % Even pair + stage(6,n+k) = stage(5,n+k) + Wn(k+1)*stage(5,n+k+16); + % Odd par + stage(6,n+k+16) = stage(5,n+k) - Wn(k+1)*stage(5,n+k+16); + end +end + +% calculations using separate real and imaginary numbers +for n = 1 : 2^5 : fftLength + for k = 0 : 15 + st_real_tmp(5,n+k+16) = ( Wr(k*1+1) * st_real(5,n+k+16) ) - ( Wi(k*1+1) * st_imag(5,n+k+16) ); + st_imag_tmp(5,n+k+16) = ( Wi(k*1+1) * st_real(5,n+k+16) ) + ( Wr(k*1+1) * st_imag(5,n+k+16) ); + end +end +for n = 1 : 2^5 : fftLength + for k = 0 : 15 + % Even pair + st_real(6,n+k) = st_real(5,n+k) + st_real_tmp(5,n+k+16); + st_imag(6,n+k) = st_imag(5,n+k) + st_imag_tmp(5,n+k+16); + % Odd par + st_real(6,n+k+16) = st_real(5,n+k) - st_real_tmp(5,n+k+16); + st_imag(6,n+k+16) = st_imag(5,n+k) - st_imag_tmp(5,n+k+16); + end +end +
Ploting out
slowly plot result
figure(5) +for n = 1 : bits +1 + %plot( abs( real_n(i, :) + j.*imag_n(i, :) ) ); + stem( abs( stage(n,1:fftLength/2) ) ); +% pause(1); +end +xt = xticks; % returns the current x-axis tick values as a vector +fstep = fs/fftLength; % tick of f axis in f domain +xtnew = round((xt-1)*fstep/1000, 1) ; % calculate new tick in kHz +xticklabels(xtnew) % set new tick labels +title('FFT using custom function') +ylabel('magnitude'), xlabel('kHz') + +figure(6) +for n = 1 : bits +1 + %plot( abs( real_n(i, :) + j.*imag_n(i, :) ) ); + temp = st_real + 1i * st_imag; + stem( abs( temp(n,1:fftLength/2) ) ); +% pause(1); +end +xt = xticks; % returns the current x-axis tick values as a vector +fstep = fs/fftLength; % tick of f axis in f domain +xtnew = round((xt-1)*fstep/1000, 1) ; % calculate new tick in kHz +xticklabels(xtnew) % set new tick labels +title('FFT using custom function real/imag separate') +ylabel('magnitude'), xlabel('kHz') + +figure(7) +dif2 = 100* abs(temp(bits +1,1:fftLength/2) - ft(1:fftLength/2))./abs(ft(1:fftLength/2)) ; +plot(dif2, 'blue') +title('Difference in calculations') +xt = xticks; % returns the current x-axis tick values as a vector +fstep = fs/fftLength; % tick of f axis in f domain +xtnew = round((xt-1)*fstep/1000, 1) ; % calculate new tick in kHz +xticklabels(xtnew) % set new tick labels +ylabel('percents, %'), xlabel('kHz') +
