% fm4_lab_exercise_C.m
%
% MATLAB function to read input and output data sequences
% saved from FM4 board in Intel hex format files
% and estimate power response of system using
% direct and indirect estimation (periodogram and correlogram)
% methods and cross correlation.
%
% Provided as starting point for B31SC Laboratory Exercise C.
%
% Students need to adjust parameter values used and plot
% results on graphs with correctly scaled and labelled axes.
%

function fm4_lab_exercise_C()
clear
% enter name of Intel hex format file containing
% input data sequence x(n), including .dat suffix
fname = input('enter input filename ','s');
fid = fopen(fname,'rt');
floatcount = 0;
dummy = fscanf(fid,'%c',1);
if (dummy ~= ':')
    disp('error: initial colon not found');
else
% process data from this file
finished = 0;
while (finished == 0)
    % move to next line
    while (fscanf(fid,'%c',1) ~= ':'); end
    % get number of 32-bit hex values on line
    N = hex2dec(fscanf(fid,'%c',2))/4;
    % read and discard next 6 characters
    fscanf(fid,'%c',6);
    if (N > 0)
        for i=1:N
          % read 8 character hex string and convert to IEEE float 754 single
          hexstring = fscanf(fid,'%c',8);
          reordered(1) = hexstring(7);
          reordered(2) = hexstring(8);
          reordered(3) = hexstring(5);
          reordered(4) = hexstring(6);
          reordered(5) = hexstring(3);
          reordered(6) = hexstring(4);
          reordered(7) = hexstring(1);
          reordered(8) = hexstring(2);
          xn(floatcount+1)= hexsingle2num(reordered);
          floatcount = floatcount+1;
        end
    else
    finished = 1;
    end
end
end
fclose(fid);          
Nx = floatcount

% enter name of Intel hex format file containing
% output data sequence y(n), including .dat suffix
fname = input('enter output filename ','s');
fid = fopen(fname,'rt');
floatcount = 0;
dummy = fscanf(fid,'%c',1);
if (dummy ~= ':')
    disp('error: initial colon not found');
else
%process data from this file
finished = 0;
while (finished == 0)
    % move to next line
    while (fscanf(fid,'%c',1) ~= ':'); end
    % get number of 32-bit hex values on line
    N = hex2dec(fscanf(fid,'%c',2))/4;
    % read and discard next 6 characters
    fscanf(fid,'%c',6);
    if (N > 0)
        for i=1:N
          % read 8 character hex string and convert to IEEE float 754 single
          hexstring = fscanf(fid,'%c',8);
          reordered(1) = hexstring(7);
          reordered(2) = hexstring(8);
          reordered(3) = hexstring(5);
          reordered(4) = hexstring(6);
          reordered(5) = hexstring(3);
          reordered(6) = hexstring(4);
          reordered(7) = hexstring(1);
          reordered(8) = hexstring(2);
          yn(floatcount+1)= hexsingle2num(reordered);
          floatcount = floatcount+1;
        end
    else
    finished = 1;
    end
end
end
fclose(fid);          
Ny = floatcount

if (Nx ~= Ny)
    disp('error: input and output sequences different lengths');
end

yvariance = 0;
for i=1:Ny
    yvariance = yvariance+yn(i)*yn(i)/Ny;
end
yvariance
xvariance = 0;
for i=1:Nx
    xvariance = xvariance+xn(i)*xn(i)/Nx;
end
xvariance

% cross correlate input and output sequences
% (2M+N) must equal less than xN and yN
N = 5000; % number of products summed for each Rxy(m) estimate
M = 200; % Rxy(m) is estimated for shift m where -M <= m <= M
for m = -M:M % estimate Rxy for +ve and -ve shifts
 Rxy(m+M+1) = 0.0; % initialise Rxx summation
                   % -M <= m <= M but index must run from
                   % 1 to (2*M + 1)
 for n = 1:N % sum N products x(n)*x(n+m)
   Rxy(m+M+1)= Rxy(m+M+1)+xn(n+M)*yn(n+M+m)/N;
 end
end
% plot representative sections of xn and yn, and (2M+1)
% calculated values of cross correlation estimate Rxy(m)
figure(1),
subplot(3,1,1),bar(xn(2000:2000+2*M+1)),title('sequence x(n)'),xlabel('n'),grid;
subplot(3,1,2),bar(yn(2000:2000+2*M+1)),title('sequence y(n)'),xlabel('n'),grid;
subplot(3,1,3),bar(-M:M,Rxy),title('cross correlation function R_{xy}(m)'),xlabel('m'),grid;
figure
% compute and plot estimate of power response based on
% squared magnitude of DFT of estimated impulse response
% divided by variance of x(n)
%
% choose section of cross correlation estimate containing
% estimated impulse response
% this sets number of DFT points
dummy = Rxy(100:360);
plot(20*log10((abs(fft(dummy)).*abs(fft(dummy)))/(xvariance*xvariance)))
figure
% compute direct estimate of power response as
% square of magnitude of NUM-point DFT of section of 
% output sequence y(n), divided by NUM times variance of x(n)
% 
NUM = 512;
dummy=yn(2000:(2000+NUM-1));
plot(20*log10((abs(fft(dummy)).*abs(fft(dummy)))/(xvariance*NUM)))
% compute indirect estimate (correlogram) of power response as
% magnitude of DFT of Bartlett-windowed section of ACF 
% of output sequence y(n), divided by variance of x(n)
% 

N = 5000; % number of products summed for each Ryy estimate
M = 200; % Ryy is estimated for shift m where -M <= m <= M
for m = -M:M % estimate Ryy for +ve and -ve shifts
 Ryy(m+M+1) = 0.0; % initialise Ryy summation
                   % -M <= m <= M but index must run from
                   % 1 to (2*M + 1)
 for n = 1:N % sum N products x(n)*x(n+m)
   Ryy(m+M+1)= Ryy(m+M+1)+yn(n+M)*yn(n+M+m)/N;
 end
end
w = bartlett(2*M+1);
Ryyw=Ryy.*w';
hw=20*log10(( abs(fft(Ryyw))/(xvariance) ));
hwf=0:2*pi/(2*M):4;

figure
plot(hwf(1:M),hw(1:M),'g')