%% Create IFFT for HDL Code Generation
%%  
% Create the specifications and input signal. This example uses a
% 128-point FFT.  
N = 128; 
Fs = 40; 
t = (0:N-1)'/Fs;
x = sin(2*pi*15*t) + 0.75*cos(2*pi*10*t);
y = x + .25*randn(size(x));
y_fixed = sfi(y,32,16);
noOp = zeros(1,'like',y_fixed);
%% 
% Compute the FFT of the signal to use as the input to the IFFT object.  
% 
% *Note:* This object syntax runs only in R2016b or later. If you are using an 
% earlier release, replace each call of an object with the equivalent |step|
% syntax. For example, replace |myObject(x)| with |step(myObject,x)|.
hdlfft = dsp.HDLFFT('FFTLength',N,'BitReversedOutput',false);
Yf = zeros(1,4*N);
validOut = false(1,4*N);
for loop = 1:1:N 
   	  [Yf(loop),validOut(loop)] = hdlfft(complex(y_fixed(loop)),true); 
end
for loop = N+1:1:4*N 
	  [Yf(loop),validOut(loop)] = hdlfft(complex(noOp),false); 
end
Yf = Yf(validOut == 1);
%% 
% Plot the single-sided amplitude spectrum.
plot(Fs/2*linspace(0,1,N/2),2*abs(Yf(1:N/2)/N))
title('Single-Sided Amplitude Spectrum of Noisy Signal y(t)')
xlabel('Frequency (Hz)')
ylabel('Output of FFT (f)')
%% 
% Select frequencies that hold the majority of the energy in the signal. The
% |cumsum| function does not accept fixed-point arguments, so convert
% the data back to |double|.
[Ysort,i] = sort(abs(double(transpose(Yf(1:N)))),1,'descend');
Ysort_d = double(Ysort);
CumEnergy = sqrt(cumsum(Ysort_d.^2))/norm(Ysort_d);
j = find(CumEnergy > 0.9, 1);
     disp(['Number of FFT coefficients that represent 90% of the ', ...
     'total energy in the sequence: ', num2str(j)])  
Yin = zeros(N,1);
Yin(i(1:j)) = Yf(i(1:j));
%%  
% Write a function that creates and calls the IFFT System object(TM). 
% You can generate HDL from this function. 
%
% <include>HDLIFFT128.m</include> 
%
%%
% Compute the IFFT by calling the function for each data sample.
Xt = zeros(1,3*N);
validOut = false(1,3*N);
for loop = 1:1:N 
    [Xt(loop),validOut(loop)] = HDLIFFT128(complex(Yin(loop)),true); 
end
for loop = N+1:1:3*N
    [Xt(loop),validOut(loop)] = HDLIFFT128(complex(0),false); 
end
%% 
% Discard invalid output samples. Then inspect the output and compare it
% with the input signal. The original input is in green.
Xt = Xt(validOut==1);
Xt = bitrevorder(Xt);
norm(x-transpose(Xt(1:N)))
figure
stem(real(Xt))
figure
stem(real(x),'--g')

%% 
% Copyright 2012 The MathWorks, Inc.