% EE 438 Digital Signal Processing with Applications 
% DTFT Application
% file: analyzer.m
% Windowing a Cosine Signal
%---------------------------------------------------------------
% Initialize
clear all;
close all;
%---------------------------------------------------------------
% Parameters
gain = 0.5;
gain = 2;	% amplitude of input signal
T = 1.0;	% duration of input, in seconds
Fs = 10000;	% sampling frequency 
Nplot = 200;	% number of samples in plots
%---------------------------------------------------------------
% Generate input signal x and put it through system
N = Fs*T;
log2N = log(N)/log(2);
ceillog2N = ceil(log2N);
N = 2^ceillog2N		% making N a power of 2
x = gain*sig_3_gen(T, Fs, N);
y = system_3(x);
%--------------------------------------------------------------
% Send signals through sound card.  Keep values between [-1,1], while
% maintaining relative amplitude.
disp('Input signal x followed by output signal y')
pause
%soundsc(resample(x,8192,Fs));
sound(x/gain,Fs);
pause
%soundsc(resample(y,8192,Fs));
sound(y/gain,Fs);
%--------------------------------------------------------------
% Generate plots
%
figure(1)
subplot(2, 1, 1), plot(x(1:Nplot))
title('System input')
xlabel('n')
ylabel('x[n]')
subplot(2, 1, 2), plot(y(1:Nplot))
title('System output')
xlabel('n')
ylabel('y[n]')
pause
%--------------------------------------------------------------
% Compute spectra
%
X = dtft(x, N);
Y = dtft(y, N);
%--------------------------------------------------------------
% Plot spectra
%
omega = -pi:2*pi/(N-1):pi;
figure(2)
subplot(2, 1, 1), plot(omega, abs(Y))
title('Output spectrum magnitude')
ylabel('|Y(\omega)|')
axis([-pi, pi, 0, max(abs(Y))])
subplot(2, 1, 2), plot(omega, abs(X))
title('Input spectrum magnitude')
ylabel('|X(\omega)|')
xlabel('\omega')
axis([-pi, pi, 0, max(abs(X))])
pause
%--------------------------------------------------------------
% SNR
%
x_sq = x.*x;
P_x = sum(x_sq)/N
e = y-x;
e_sq = e.*e;
P_e = sum(e_sq)/N
SNR = 10*log10(P_x/P_e)
%--------------------------------------------------------------
% Error spectrum
%
E = dtft(e, N);
figure(3)
subplot(2, 1, 1), plot(e(1:Nplot))
title('Error')
xlabel('n')
ylabel('e[n]')
omega = -pi:2*pi/(N-1):pi;
subplot(2, 1, 2), plot(omega, abs(E))
title('Error spectrum magnitude')
ylabel('|E(\omega)|')
xlabel('\omega')
if max(abs(E))~=0
  axis([-pi, pi, 0, max(abs(E))])
end
pause
%--------------------------------------------------------------
% Null error at fundamental frequencies
%
E_max = max(abs(E));
epsilon = 0.01*E_max;
i_max = find(abs(E) >= E_max - epsilon);
delta = i_max(2) - i_max(1);
mask = ones(1,N);
null = zeros(1,2*delta);
null_start = N/2-delta;
null_stop = null_start + length(null) - 1;
mask(null_start:null_stop) = null;
E_h = mask.*E;
figure(4)
plot(omega, abs(E_h))
title('Error harmonics only')
ylabel('|E_h(\omega)|')
xlabel('\omega')
if max(abs(E_h))~=0
  axis([-pi, pi, 0, max(abs(E_h))])
end 
X = dtft(x,N);
%--------------------------------------------------------------
% Total harmonic distortion (THD)
%
X = dtft(x,N);
X_sq = abs(X).^2;
P_X = sum(X_sq)/N
E_sq = abs(E_h).^2;
P_E_h = sum(E_sq)/N
total_harmonic_distortion = (P_E_h/P_X)*100