% -------------------------------------------------------------------------
% run_g7play_blockbyblock.m (edited 1/13/16)

% Matlab program to run the MU ("modified unit") box which simultaneously 
% outputs voltages on 2 channels, and measures voltage input on 2 channels 

% IMPORTANT NOTE: This program is designed to run from inside a folder at
% the same level as data, bin, WINDOWS, LINUX, MAC, etc. All of the calls
% to files are programs are written with this in mind. You must change
% these calls if you move the location of this file.
% -------------------------------------------------------------------------

addpath ../bin;

% Choose one of the following (g7string plays the executable from the
% correct directory). NOTE: Test the command "system(g7string)" in the 
% command line from the same folder as this function 
% (g7play should run if the string is correct).

%addpath MAC; g7string = '../MAC/g7play';
addpath WINDOWS; g7string = '..\WINDOWS\g7play';
%addpath LINUX; g7string = '../LINUX/g7play';
 
% -------------------------------------------------------------------------
% Initialize parameters & settings
% -------------------------------------------------------------------------
NB = 100; % number of blocks
separate_blocks = 1; 
%0: save & plot average of all blocks
%1: save and plot individual blocks --- 'block by block mode'!!!
% in 'block by block' mode you can look at the data in quasi-real-time (lags
% far behind real time)

NFT=1024; % discrete fourier transform length 
NF=NFT/2+1; % number of frequency bins
Fs=48000; % sampling frequency
Ts=1/Fs; % sampling period
Fmax=Fs/2; % max frequency
Fmin=Fs/NFT; % minimum (non-zero) frequency & frequency bin size
f=Fmin*(0:NF-1); % vector of frequencies (in Hertz)
t=Ts*(0:NFT-1); % vector of times (in seconds)
G=0.01; % gain paramter
% NOTE: Time-domain signal is normalized so that the maximum value is 1. 
% The gain value G scales this signal. Example: If G =0.01, then the
% maximum value of the time-domain signal will be 0.01.

global V,G; % global parameters to be read by rd1.m
V=struct( 'Fs',Fs, 'NFT',NFT, 'f',f, 't',t); 

%Basic parameters for g7play
Nmax=2^31-1; % max integer 
% Navg doesn't get used in this program, since we do block-by-block
%Navg=32; % number of blocks sent/received to be averaged together
%Navg must be a power of 2, for ola.m (overlap add) to work!

% -------------------------------------------------------------------------
% Generate output signals (left and right channel)
% -------------------------------------------------------------------------

phi=(0:NF-1)/(NF-1); % initialize vector for chirp phase
phi=(NFT/2)*pi*phi.*phi; % set frequency-domain phase of chirp (quadratic) 
Xf=exp(-1j*phi); % frequency-domain signal
xt=fsst(Xf); % calculate time-domain signal (fsst.m takes the inverse fourier transform of C(f))
xt=0.999*xt/max(abs(xt)); % normalize time-domain signal so its highest value is almost 1
Xf=fast(xt); % re-calculate frequency-domain signal (fast.m takes the fourier transform of c(t))

Lf=Xf; % set frequency-domain signal for "left" output channel
Rf=Xf; % set frequency-domain signal for "right" output channel
lt=fsst(Lf); rt=fsst(Rf); % calculate corresponding "time-domain" signals

% apply gain parameter
Gr=G; Gl=G; % set channel gain equal to G for both channels
lt=round(Gl*lt*Nmax); %Scale data to int32 2^31-1 range & apply gain
rt=round(Gr*rt*Nmax);

%save the data in files stim0 and stim1
fd=fopen('../data/stim0','w');
N=fwrite(fd,lt,'int32'); %Scale to Nmax (32 bit integer)
fclose(fd);

fd=fopen('../data/stim1','w');
N=fwrite(fd,rt,'int32'); %Scale to Nmax (32 bit integer)
fclose(fd);

delete('../data/resp*') % delete any previously stored responses

h= figure();
h1 = subplot(221);
h2 = subplot(223);
h3 = subplot(222);
h4 = subplot(224);

P = get(h,'Position'); P(1) = P(1) - P(3)/2; P(3) = 2*P(3); set(h,'Position',P);

c0 = zeros(NB,NFT);
c1 = zeros(NB,NFT);
%% Measure & Plot

% -Cx command saves all blocks to 1 file
% -Ux command saves blocks to individual files

if separate_blocks ==0
    command =sprintf(['%s -Q 0 -B ',num2str(NB),' -R0 ../data/stim0 -R1 ../data/stim1 -C0 ../data/resp0 -C1 ../data/resp1 -S %d &'],g7string,NFT); % chirp
    system(command);
  
    fd=fopen('../data/resp0'); c0=fread(fd,'int32'); fclose(fd);
    fd=fopen('../data/resp1'); c1=fread(fd,'int32'); fclose(fd);
    
    c0=ola(c0,NFT); %overlap-add scaled data
    c1=ola(c1,NFT); %overlap-add scaled data
    C0=fast(c0);
    C1=fast(c1);

    figure(h)
    %subplot(221)
    plot(h1,(0:NFT-1)*(1/Fs)*1000,c0/Nmax,'k','linewidth',1); %hold(gca,'off')
    ylabel(h1,'Volts'); title(h1,'Voltage Waveform (Response Channel 0)'); xlabel(h1,'Time [ms]')
    %subplot(223)
    plot(h2,(0:NFT-1)*(1/Fs)*1000,c1/Nmax,'k','linewidth',1); %hold(gca,'off')
    ylabel(h2,'Volts'); title(h2,'Voltage Waveform (Response Channel 1)'); xlabel(h2,'Time [ms]')
    %subplot(222)
    loglog(h3,f/10^3,abs(C0)/Nmax,'k') ;
    title(h3,'Voltage Spectrum (Response Channel 0)'); ylabel(h3,'Volts')
    xlabel(h3,'Frequency [kHz]'); set(h3,'ylim',[10^-4 10^2]);
    %subplot(224)
    loglog(h4,f/10^3,abs(C1)/Nmax,'k') ;
    title(h4,'Voltage Spectrum (Response Channel 1)'); ylabel(h4,'Volts')
    xlabel(h4,'Frequency [kHz]'); set(h4,'ylim',[10^-4 10^2]);
else
    command =sprintf(['%s -Q 0 -B ',num2str(NB),' -R0 ../data/stim0 -R1 ../data/stim1 -U0 ../data/resp0 -U1 ../data/resp1 -S %d'],g7string,NFT); % chirp
    system(command);
   
    for k1 = 1:NB
        
        %turn this on if the plot gets ahead of g7play (this probably won't
        %happen)
        %while flag ~=2, flag = exist(['data/resp1.',num2str(k1)],'file'); end % don't get ahead of ourselves
        %flag = 0;
                        
        fd=fopen(['../data/resp0.',num2str(k1)]); c0(k1,:)=fread(fd,'int32'); fclose(fd);
        fd=fopen(['../data/resp1.',num2str(k1)]); c1(k1,:)=fread(fd,'int32'); fclose(fd);
        
        figure(h)
        %subplot(221)
        plot(h1,(0:NFT-1)*(1/Fs)*1000,c0(k1,:)/Nmax,'k','linewidth',1); %hold(gca,'off')
        ylabel(h1,'Volts'); title(h1,'Voltage Waveform (Response Channel 0)'); xlabel(h1,'Time [ms]')
        %subplot(222)        
        plot(h2,(0:NFT-1)*(1/Fs)*1000,c1(k1,:)/Nmax,'k','linewidth',1); %hold(gca,'off')
        ylabel(h2,'Volts'); title(h2,'Voltage Waveform (Response Channel 1)'); xlabel(h2,'Time [ms]')
        
        n=NFT/2+1;
        C0=fft(c0(k1,:));
        C0=C0(1:fix(n));
        C1=fft(c1(k1,:));
        C1=C1(1:fix(n));
        
        %subplot(223)
        loglog(h3,f/10^3,abs(C0)/Nmax,'k') ;
        title(h3,'Voltage Spectrum (Response Channel 0)'); ylabel(h3,'Volts')
        xlabel(h3,'Frequency [kHz]'); set(h3,'ylim',[10^-4 10^2]);
        %hold(gca,'off')

        loglog(h4,f/10^3,abs(C1)/Nmax,'k') ;
        title(h4,'Voltage Spectrum (Response Channel 1)'); ylabel(h4,'Volts')
        xlabel(h4,'Frequency [kHz]'); set(h4,'ylim',[10^-4 10^2]);

    end
end