Create an i.i.d White Gaussian Vector X

% create x
Num = 1000;
x = randn(1,Num);

PMF, Autocorrelation, and PSD of X

% plot the normalized histogram of x
[count,center] = hist(x,50); %50 bins
figure; stem(center,count./Num); axis tight;
title('PMF of X')

% plot the autocorrelation of x
Rx = xcorr(x);
figure; plot([-Num+1:1:Num-1],Rx); axis tight;
title('Autocorrelation of X')

% plot the spectral density of x
Lfft = 1024; % next power of 2 from length of signal
Sx = fftshift(abs(fft(Rx,Lfft)));
figure; plot([-pi :2*pi/(Lfft-1): pi],Sx); axis tight;
title('Power Spectral Density of X')   Impulse and Frequency Response of H

% impulse response h[n]
impulse = zeros(1,10);
impulse(1) = 1;
h = filter(1,[1 -1/3],impulse);
figure; stem(h); axis tight;
title('Impulse Response h[n]')

% frequency response H(\Omega) using freqz
figure; freqz(1,[1 -1/3]); axis tight;
title('Frequency Response H(\Omega)')

% magnitude of frequency response using fft of h[n]
figure; plot([-pi :2*pi/(Lfft-1): pi], fftshift(abs(fft(h,Lfft)))); axis tight;
title('Frequency Response H(\Omega)')   X -> H -> Y

y = filter(1,[1 -1/3],x);
figure; plot(abs(x-y)); axis tight;
title('Absolute Difference Between X and Y') PMF, Autocorrelation, and PSD of Y

% plot the normalized histogram of y
[count,center] = hist(y,50); % 50 bins
figure; stem(center,count./Num); axis tight;
title('PMF of Y');

% plot the autocorrelation of x
Ry = xcorr(y);
figure; plot([-Num+1:1:Num-1],Ry); axis tight;
title('Autocorrelation of Y')

% plot the spectral density of x
Lfft = 1024; % next power of 2 from length of signal
Sy = fftshift(abs(fft(Ry,Lfft)));
figure; plot([-pi :2*pi/(Lfft-1): pi],Sy); axis tight;
title('Power Spectral Density of Y')   Y -> Hinv -> Xhat

xhat = filter([1 -1/3],1,y); %reverse coefficients for hinv
figure; plot(abs(x-xhat)); axis tight;
title('Absolute Difference Between X and Xhat')

% frequency response H(\Omega) using freqz
figure; freqz([1 -1/3],1); axis tight;
title('Frequency Response H^{-1}(\Omega)')

% autocorrelation/psd of xhat (should be white again)
Rxhat = xcorr(x);
Lfft = 1024; % next power of 2 from length of signal
Sxhat = fftshift(abs(fft(Rxhat,Lfft)));

% compare spectral densities
figure;

subplot(1,3,1); plot([-pi :2*pi/(Lfft-1): pi],Sx); axis tight;
title('Power Spectral Density of X');

subplot(1,3,2); plot([-pi :2*pi/(Lfft-1): pi],Sy); axis tight;
title('Power Spectral Density of Y');

subplot(1,3,3); plot([-pi :2*pi/(Lfft-1): pi],Sxhat); axis tight;
title('Power Spectral Density of \hat{X}');   