% -------------------------------------------------------------------------
% show_data.m (created 1/13/17)
%
% Program to read in the data written by run_g7play
% ------------------------------------------------------------------------- 
% NOTE: In general, upper case signal variables represent frequency domain,
% while lowercase signal variables represent time domain

% FLAGS: 0 = no, 1 = yes (boolean)
show_stimulus = 0; %FLAG to plot stimulus
plot_impedance = 1; % FLAG to plot impedance (only relevant for a particular
% circuit - see comments and lab manual)
dechirp = 0; % FLAG (make this zero if your stimulus wasn't a chirp)

global V;

% close all % close open figs (NOTE: If too many figs are open, Matlab will
% crash)

try % running from run_g7play.m in directory above
    load data/V.mat
    fd=fopen('data/stim0'); s0=fread(fd,'int32'); fclose(fd);
    fd=fopen('data/stim1'); s1=fread(fd,'int32'); fclose(fd);
    fd=fopen('data/resp0'); c0=fread(fd,'int32'); fclose(fd);
    fd=fopen('data/resp1'); c1=fread(fd,'int32'); fclose(fd);
catch % running from this directory
    load ../data/V.mat
    fd=fopen('../data/stim0'); s0=fread(fd,'int32'); fclose(fd);
    fd=fopen('../data/stim1'); s1=fread(fd,'int32'); fclose(fd);
    fd=fopen('../data/resp0'); c0=fread(fd,'int32'); fclose(fd);
    fd=fopen('../data/resp1'); c1=fread(fd,'int32'); fclose(fd);
end

Nmax=2^31-1;

% drop the first block and calculate response in Volts
c0 = V.VmaxAD*c0(V.NFT+1:end)/Nmax;
c1 = V.VmaxAD*c1(V.NFT+1:end)/Nmax;
% average the blocks
c0=mean(reshape(c0(:),V.NFT,length(c0(:))/V.NFT),2);
c1=mean(reshape(c1(:),V.NFT,length(c1(:))/V.NFT),2);
% calculate frequency-domain responses
C0=fast(c0);
C1=fast(c1);

% calculate stimulus in Volts
s0 = V.VmaxDA*s0/Nmax; S0 = fast(s0);
s1 = V.VmaxDA*s1/Nmax; S1 = fast(s1);

% dechirp responses if flag indicates the stimulus was a chirp
if dechirp
    C0 = C0.*exp(-1j*unwrap(angle(S0))); c0 = fsst(C0);
    C1 = C1.*exp(-1j*unwrap(angle(S1))); c1 = fsst(C1);
end

% initialize figures
h1= figure(); set(h1,'units','normalized','outerposition',[.03 .15 .9 .7]);

% Plot time and frequency responses from both channels
figure(h1)
subplot(231)
plot(V.t*1000,c0,'k','linewidth',1);
ylabel('Volts'); title('Voltage Waveform (RESPONSE Channel 0)'); xlabel('Time [ms]')
subplot(234)
plot(V.t*1000,c1,'k','linewidth',1);
ylabel('Volts'); title('Voltage Waveform (RESPONSE Channel 1)'); xlabel('Time [ms]')
subplot(232)
loglog(V.f/10^3,abs(C0),'k') ;
title('Voltage Spectrum Magnitude (RESPONSE Channel 0)'); ylabel('Volts')
xlabel('Frequency [kHz]'); set(gca,'ylim',[10^-4 10^2]);
subplot(235)
loglog(V.f/10^3,abs(C1),'k') ;
title('Voltage Spectrum Magnitude (RESPONSE Channel 1)'); ylabel('Volts')
xlabel('Frequency [kHz]'); set(gca,'ylim',[10^-4 10^2]);
subplot(233)
semilogx(V.f/10^3,unwrap(angle(C0))/pi,'k') ;
title('Voltage Spectrum Phase (RESPONSE Channel 0)'); ylabel('Radians/\pi')
xlabel('Frequency [kHz]');
subplot(236)
semilogx(V.f/10^3,unwrap(angle(C1))/pi,'k') ;
title('Voltage Spectrum Phase (RESPONSE Channel 1)'); ylabel('Radians/\pi')
xlabel('Frequency [kHz]');

if plot_impedance == 1 % Plot the impedance
    % Voltage divider: Z_{test} + R_{ref} in series.
    % Z_{test} is a unknown impedance - resistor, capacitor, loudspeaker, etc.
    % For the code below (C0=V1, C1=V2), where we are assuming V1 and V2 are
    % measuring across Z_{test} (V2 is between Z_{test} and R_{ref}).
    % NOTE: Make sure your answer makes sense - the channels may easily be swapped by accident!
    
    Rref=22; % voltage divider resistor
    Zx=Rref*(C0./C1 -1);
    
    h2 = figure();
    subplot(211)
    loglog(V.f/1000,abs(Zx),'b', 'linewidth',2);
    title(['Impedance Magnitude [Z = Rref*(V1/V2-1)], R_{ref} = ',num2str(Rref),' \Omega']);
    xlabel('Frequency [kHz]'); ylabel('|Z| [Ohm]')
    grid on
    subplot(212)
    semilogx(V.f/1000,unwrap(angle(Zx))/pi,'b', 'linewidth',2);
    title(['Impedance Phase [Z = Rref*(V1/V2-1)], R_{ref} = ',num2str(Rref),' \Omega']);
    ylabel('\angle Z [Rad/\pi]')
    xlabel('Frequency [kHz]');
    grid on
end

if show_stimulus == 1
    h3= figure(); set(h3,'units','normalized','outerposition',[.08 .05 .9 .7]);
    
    % Plot time and frequency stimuli from both channels
    figure(h3)
    subplot(231)
    plot(V.t*1000,s0,'k','linewidth',1);
    ylabel('Volts'); title('Voltage Waveform (STIMULUS Channel 0)'); xlabel('Time [ms]')
    subplot(234)
    plot(V.t*1000,s1,'k','linewidth',1);
    ylabel('Volts'); title('Voltage Waveform (STIMULUS Channel 1)'); xlabel('Time [ms]')
    subplot(232)
    loglog(V.f/10^3,abs(S0),'k') ;
    title('Voltage Spectrum Magnitude (STIMULUS Channel 0)'); ylabel('Volts')
    xlabel('Frequency [kHz]'); set(gca,'ylim',[10^-4 10^2]);
    subplot(235)
    loglog(V.f/10^3,abs(S1),'k') ;
    title('Voltage Spectrum Magnitude (STIMULUS Channel 1)'); ylabel('Volts')
    xlabel('Frequency [kHz]'); set(gca,'ylim',[10^-4 10^2]);
    subplot(233)
    semilogx(V.f/10^3,unwrap(angle(S0))/pi,'k') ;
    title('Voltage Spectrum Phase (STIMULUS Channel 0)'); ylabel('Radians/\pi')
    xlabel('Frequency [kHz]');
    subplot(236)
    semilogx(V.f/10^3,unwrap(angle(S1))/pi,'k') ;
    title('Voltage Spectrum Phase (STIMULUS Channel 1)'); ylabel('Radians/\pi')
    xlabel('Frequency [kHz]');
end