function myFFTTest(N1,sRt,dur,sigProb,fF, thinningProb)
%
% example: myFFTTest(6*1024,44000,2,0.0005,30,0.3)

%---------------------------------------------------------%
% CONSTANTS
%---------------------------------------------------------%
fMin =0;    %Minimum in shown Spectrum in Hz
fMax = 200; %Maximum in shown Spectrum in Hz
%sigProb=0.005; percentage of signal-occurences in the random signal
%N1=4*2048; % Number of FFT points
%sRt = 42000 % sampling Rate in Hz
%dur=3 %duration of the signal in sec
%fF=30; %particular Frequency for non-random signal generation

    
close all;

%---------------------------------------------------------%
% generate time Vector
%---------------------------------------------------------%
t=(0:1/sRt:dur);               
%---------------------------------------------------------%
% generate spike-signal
%---------------------------------------------------------%
xSig=zeros(1,length(t)); %generate a zero vector
xSig(rand(1,length(t))<=sigProb)=1; %set random indices to 1




%---------------------------------------------------------%
% FAST FOURIER TRANSFORM of the original signal
%---------------------------------------------------------%
f=(0:N1-1)*sRt/N1;   %Defining the frequency points [axis]
XR=abs(fft(xSig,N1));  %find the magnitude of the FFT using No windoing
    %Windowing: Hann
        xh1=hamming(size(xSig,2)); %Define the hamming samples
        xw=xSig.*xh1';  % Window the signal
        XH=abs(fft(xw,N1));  %find the magnitude of the FFT of the windowed signal
    %Windowing Hamming
        xh2=hann(size(xSig,2)); %Define the HANN samples
        xw2=xSig.*xh2';  % Window the signal
        XH2=abs(fft(xw2,N1));  %find the magnitude of the FFT of the windowed signal
        
        %---------------------------------------------------------%
        % Graphical Output
        %---------------------------------------------------------%
        strTitle=sprintf('Random Noise with FFT; %i occurences',sum(xSig));
        figure('Name',strTitle,'NumberTitle','on');
            %Plotting the signal
            subplot(2,1,1); %<---??? but I need the handle for setting the xLim
                plot(t,xSig);
                title('Random signal x(t)');
                grid;
                xlabel('Time, sec');
                ylabel('Signal');

            %Plotting the FFT
            h5=subplot(2,1,2); %<---??? but I need the handle for setting the xLim
                plot(f(1:N1/2),20*log10(XR(1:N1/2)/max(XR)),f(1:N1/2),20*log10(XH(1:N1/2)/max(XH)),f(1:N1/2),20*log10(XH2(1:N1/2)/max(XH2)));
                set(h5,'xLim',[fMin fMax]); % show only a small area of the full spectrum
                title('FFT Spectrum of random signal x(t) (Rectangular, Hann, and Hamming Windows)');
                grid;
                legend('Rectangular','Hamming','Hann');
                xlabel('Frequency, Hz');
                ylabel('Normalised Magnitude, [dB]');


                
%---------------------------------------------------------%
% FAST FOURIER TRANSFORM of a particular signal
%---------------------------------------------------------%
    xSig2=zeros(1,length(t)); %create zero-vector
    xSig2(1:round(sRt/fF):end) = 1;  %fill with particular signal
    %thin out the signal
    iO=sum(xSig2);
    ctr=0;
    for i=1:sum(xSig2)
        if rand<=thinningProb  %thin out signal according to chance
            if (i*round(sRt/fF))+1 <=length(xSig2) %ensure that the vector is not enlongated
                xSig2((i*round(sRt/fF))+1)=0;
                ctr=ctr+1;
            end
        end
    end
    strTitle=sprintf('Particular signal; %i occurences (= %i %% persistency)',iO-ctr,100-round((100*ctr)/iO));
    figure('Name',strTitle,'NumberTitle','on');
    [xFrq y1 y2 y3]=myFFT(t,xSig2,N1, sRt, fMin, fMax);
    
    
%---------------------------------------------------------%
% FAST FOURIER TRANSFORM of a particular signal + random noise
%---------------------------------------------------------%
    xSig2(xSig==1)=1;        %add random noise to particular signal
    figure('Name','Particular signal + random noise with FFT','NumberTitle','on');
    [xFrq yy1 yy2 yy3]=myFFT(t,xSig2,N1, sRt, fMin, fMax);
    
    figure('Name','Difference between "noise+signal" and "noise only"','NumberTitle','on');
    %plot(xFrq,yy1-y1,xFrq,yy2-y2,xFrq,yy3-y3);
    plot(f(1:N1/2),yy1-20*log10(XR(1:N1/2)/max(XR)),f(1:N1/2),yy2-20*log10(XH(1:N1/2)/max(XH)),f(1:N1/2),yy3-20*log10(XH2(1:N1/2)/max(XH2)));
        xlim([fMin fMax]); %ylim([-100 0]);
        title('Signal Difference made by adding a particular frequency to the noise(Rectangular, Hann, and Hamming Windows)');
        grid;
        legend('Rectangular','Hamming','Hann');
        xlabel('Frequency, Hz');
        ylabel('Normalised Magnitude, [dB]');
end


%---------------------------------------------------------%
% HELPING FUNCTIONS: keeping the rest of the code simpler
%---------------------------------------------------------%

function [xFrq y1 y2 y3]=myFFT(t,xSig,N1, sRt, fMin, fMax)

% xSig: signal to be analysed
% N1: number of FFT points
% sRt: sampling rate for signal detection
% strTitle: string variable containing the title
% fMin: minimum value shown on x-Axis
% fMax: maximum value shown on x-Axis
% plotAverage: adds an optional FFT with all three windows averaged
    
    %----------------------------------------------------%
    % Calculation steps
    %----------------------------------------------------%
        f=(0:N1-1)*sRt/N1;   %Defining the frequency points [axis]
        %Windowing: None
            XR=abs(fft(xSig,N1));  %find the magnitude of the FFT using No windoing
        %Windowing: Hann
            xh1=hamming(size(xSig,2)); %Define the hamming samples
            xw1=xSig.*xh1';  % Window the signal
            XH1=abs(fft(xw1,N1));  %find the magnitude of the FFT of the windowed signal
        %Windowing Hamming
            xh2=hann(size(xSig,2)); %Define the HANN samples
            xw2=xSig.*xh2';  % Window the signal
            XH2=abs(fft(xw2,N1));  %find the magnitude of the FFT of the windowed signal
        xFrq=f(1:N1/2);
        y1=20*log10(XR(1:N1/2)/max(XR));
        y2=20*log10(XH1(1:N1/2)/max(XH1));
        y3=20*log10(XH2(1:N1/2)/max(XH2));
        
        %---------------------------------------------------------%
        % Graphical Output
        %---------------------------------------------------------%
            %Plotting the signal
            subplot(2,1,1);
                plot(t,xSig);
                title('Signal x(t)');
                grid;
                xlabel('Time, sec');
                ylabel('Signal');
            %Plotting the FFT
            subplot(2,1,2);
                plot(xFrq,y1,xFrq,y2,xFrq,y3);
                xlim([fMin fMax]); %ylim([-100 0]);
                title('FFT Spectrum of Signal (Rectangular, Hann, and Hamming Windows)');
                grid;
                legend('Rectangular','Hamming','Hann');
                xlabel('Frequency, Hz');
                ylabel('Normalised Magnitude, [dB]');
end


            
function    out=makeEmFitFun(src,fF,altFP)
% This function artificially fits the data in "src" to the given and
% weighted 2-D-Frequency-matrix "fF":
% "fF" example:   fF=[30 50 90; 5 2 1];
% "altFP": [0..1]: denoting the probability for changing a value


%--------------------------------------------------------------------% 
 %ensure that 2nd row is relative weight in the range of [0..1]
 %--------------------------------------------------------------------% 
 fF(2,:)= fF(2,:)/sum( fF(2,:));
 
 %--------------------------------------------------------------------% 
 %fill a vector according to weigthed frequencies:
 %--------------------------------------------------------------------% 
 uMax=30;
 for k=1:length(fF)
     if k==1
         x=repmat(fF(1,k),round(uMax*fF(2,k)),1)';
     else
         x=[x repmat(fF(1,k),round(uMax*fF(2,k)),1)'];
     end
 end
%--------------------------------------------------------------------% 
%correct the real values, such that they fit to the frequencies in a
%weighted manner
%--------------------------------------------------------------------% 
srcD=zeros(size(src)); %
rE=zeros(1,numel(src));
 for i=1:numel(src)
     if ~isnan(src(i))
         vN=length(fF); %amount of the frequencies to be fitted
         x=x(1,randperm(length(x)));  % "shaking" the box
         F=1/x(round((rand*(length(x)-1)))+1);     % "pulling" the lottery cards; the considered frequencies
         rE(i)=1/F;
         d=mod(src(i),F); %differencial matrix containing the correction values; for inspection purpose only
        if src(i)<F
            srcD(i)=F-src(i);

        elseif src(i)>F
            if d>F/2
                srcD(i)=F-d;
            else
                srcD(i)=-d;
            end
        end
        %---------------------------------------------------%
        % data is changed only with a certain probability
        %---------------------------------------------------%        
        if rand>altFP
            src(i)=src(i)+srcD(i);
        else
            rE(i)=-1; % if it's not changed, this is counted
        end
     end
 end
nonAltNum=abs(sum(rE(rE==-1))); %counting the non-altered entries (thus being original still)
%totNum=numel(src(isnan(src)));
totNum=numel(src);
altFrac=1-(nonAltNum/totNum); %altered fraction; facts
rE(rE==-1)=[]; % delete the non-changed
rE(rE==0)=[]; % delete the empty places
%sum(srcD)
 'Warning: results altered; output for inspection purposes only'
 
%--------------------------------------------------------------------% 
%Graphical output
%--------------------------------------------------------------------% 
    hold off;
    strTitle = sprintf('Altered Results %d out of %d (%d percent)',totNum-nonAltNum, totNum,altFrac*100);
    %strTitle = sprintf('%i percent artificially altered results: used frequencies',altFrac);
    fig2=figure('Name',strTitle,'NumberTitle','on');
    scrsz = get(0,'ScreenSize'); % get screensize
    set(fig2,'OuterPosition',[20 scrsz(4)/4.5 scrsz(3)/2 scrsz(4)/2]); %LinkerRand UntererRand Weite Höhe

    if length(rE)==1
        hist(rE); %    hist(rE,length(unique(rE)));
    elseif length(rE)>1
        hist(rE,unique(rE)); %    hist(rE,length(unique(rE)));
    end
    title(strTitle);
    xlabel('Frequencies in Hz');
    ylabel('How often a frequency was considered');
%--------------------------------------------------------------------% 
% Return Value
%--------------------------------------------------------------------% 
    out=src;
end