% This function generates a two-sided FFT plot including some aliases and
% returns the figure handle
%
%   Author: Ivo Grondman
%
function figh = fftplot(y,Fs,N,Fs2)
% If argument Fs2 is missing, set to Fs
if ~exist('Fs2','var')
    Fs2 = Fs;
end

% Fetch the length of the signal
L = length(y);

% Take the absolute value of the fourier transform (scaling by L is
% required in MATLAB)
FFTy2 = abs(fft(y)/L);

% Calculate the frequencies
df = Fs/L;
f = (0:L-1)*df;

% The frequencies larger than Fs/2 should be reduced by Fs
f(f > Fs/2) = f(f > Fs/2) - Fs;

% Shift axes 
f = fftshift(f);
FFTy2 = fftshift(FFTy2);

% Create aliases
fleft = [];
fright = [];
for k = 1:N
    fleft = [f-k*Fs2 fleft];
    fright = [fright (f+k*Fs2)];
end
falias = repmat(FFTy2,N,1);

figh = plot(f,FFTy2,fleft,falias,'r',fright,falias,'r');
xlabel('Hz');