%% cscn_epr_filtering.m.m
% This script illustrates that when filtering a CSCN signal by nulling 
% (large) portions of the FFT, the signal remains CSCN.
% Furthermore the limited steepnes of the analysis window when performing 
% PSD estimation is illustrated: In the FFT we see the hard nulling,
% whereas in the PSD there is a smoother tranisition visible.

% Tested with Matlab R2015a

% Karl Freiberger, 25. Nov. 2016
% Graz University of Technology, SPSC Lab

%% Clean up
clc;
close all;
clear all;

%% Generate Signals
nS = 128e3; % number of samples 
% 128e3 = 320(symbols)*4e-6 (symbolduration) at 100 MHz rate.
vS = randn(nS,1) + 1i*randn(nS,1); % Complex Gauss source signal

vIndexNull = (1:nS).';
nBands = 128; % number of equal bands

mIndexNull = reshape(vIndexNull, nS/nBands, nBands);
vBandsNotNulled = [10, 50, 100, 102].'; % all but these bands are nulled
% vBandsNotNulled = [10:10:nBands].'; % all but these bands are nulled
% vBandsNotNulled = [1:64, 66:nBands].'; % only band 65 is nulled
sPercentNulled = (1-length(vBandsNotNulled)/nBands)*100;
disp(sprintf('%2.2f Percent of the bandwidth are nulled', sPercentNulled))

mIndexNull(:,vBandsNotNulled) = []; 
vIndexNull = mIndexNull(:);

vSf = fft(vS);
vXf = vSf;
vXf(vIndexNull) = 0;
vX = ifft(vXf);

vX = vX./sqrt(mean(abs(vX).^2)).*sqrt(mean(abs(vS).^2)); % gain scaling

%% FFT Plot 
sEpsi = 1e-5; % regularization --> 1e-5 leads to -100 dB lower floor
figure;
  plot(20*log10(abs(vSf)+sEpsi));
  hold on;
  plot(20*log10(abs(vXf)+sEpsi));
  xlabel('FFT Frequency Bin');
  legend('Input signal FFT', ...
    'Output signal FFT = input Signal FFT with large parts nulled');
 
%% PSD Plot 
nWin = 1024*20;
hWin = @hann;
figure;
  pwelch(vS, hWin(nWin), nWin/2);  
  hold all;
  pwelch(vX, hWin(nWin), nWin/2); 

%% Time-Domain Plot 
figure;
  subplot(411);
    plot(real([vS,vX]));
    title('Real Part');
    legend('Input s[n]', 'Output x[n]')
  subplot(412);
    plot(imag([vS,vX]));
    title('Imaginary Part');
  subplot(413);
    plot(abs([vS,vX]));
    title('Absolute Value');
  subplot(414);
    plot(angle([vS,vX])./pi);
    title('Angle');
    ylabel('rad/ \pi')
  xlabel('Discrete Time Index n');
  
%% CCDF Plot
vVal = linspace(min(abs([vS;vX])), max(abs([vS;vX])),512).';  
vData = abs(vS);
vCdfS = histc(-vData,[-inf;-flipud(vVal)]);
vCdfS = flipud(cumsum(vCdfS(1:end-1)))/nS;    

vData = abs(vX);
vCdfX = histc(-vData,[-inf;-flipud(vVal)]);
vCdfX = flipud(cumsum(vCdfX(1:end-1)))/nS;    

figure;
  semilogy(vVal, vCdfS); hold all;
  semilogy(vVal, vCdfX);
  ylabel('CCDF');
  xlabel('Absolute value');

%% PDF Plot
vVal = linspace(min(abs([vS;vX])), max(abs([vS;vX])),128).'; 
hPdf = @(vData, vVal) hist(vData, vVal) / trapz(vVal, hist(vData, vVal));
nBins = 128;
figure;
  subplot(411);
    hFun = @real;
    vVal = linspace(min(hFun([vS;vX])), max(hFun([vS;vX])),nBins).';
    plot(vVal,hPdf(hFun(vS), vVal)); 
    hold all;
    plot(vVal,hPdf(hFun(vX), vVal)); 
    xlabel('Real Part');
    legend('Input s[n]', 'Output x[n]')
    title ('Histogram-Based Probabilty Density Function (PDF) Estimate ');
  subplot(412);
    hFun = @imag;
    vVal = linspace(min(hFun([vS;vX])), max(hFun([vS;vX])),nBins).';
    plot(vVal,hPdf(hFun(vS), vVal)); 
    hold all;
    plot(vVal,hPdf(hFun(vX), vVal));  
    xlabel('Imaginary Part');
  subplot(413);
    hFun = @abs;
    vVal = linspace(min(hFun([vS;vX])), max(hFun([vS;vX])),nBins).';
    plot(vVal,hPdf(hFun(vS), vVal)); 
    hold all;
    plot(vVal,hPdf(hFun(vX), vVal)); 
    xlabel('Absolute Value');
  subplot(414);
    hFun = @angle;
    vVal = linspace(min(hFun([vS;vX])), max(hFun([vS;vX])),nBins).';
    plot(vVal,hPdf(hFun(vS), vVal)); 
    hold all;
    plot(vVal,hPdf(hFun(vX), vVal)); 
    xlabel('Angle (rad)');
  

%% Hypothesis test.
% H0 is that vS and vX are from the same distribution. 
% If they are from the same distribution, kstest2 returns 0
% if the distributions kstest2 returns 1
% Alpha is the significance intervall (a low value means very high
% significance)
 [bDiffDistr, sPValue] = kstest2(abs(vS),abs(vX), 'Alpha', 0.05)
