% Matlab program to test the simultaneously measuring on 2 channels,

clear all;
close all;

addpath('bin');
%addpath ./MAC;
addpath('WINDOWS');
%addpath ./LINUX;
% if the program crashes with message "'g7play' is not recognized as an
% internal or external command" you should put g7play in the same directory
% as this program

NFT=1024; NF=NFT/2+1; Fs=48000; Ts=1/Fs;
Fmax=Fs/2; Fmin=Fs/NFT; f=Fmin*(0:NF-1); 
t=Ts*(0:NFT-1);
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)

global V; V=struct( 'Fs',Fs, 'NFT',NFT, 'f',f, 't',t); 

%Basic parameters for g7play
Nmax=2^31-1; Navg=4; %This must be a power of 2, for ola() to work!

%Generate chirp for left and right channels
fee=(0:NF-1)/(NF-1);
fee=(NFT/2)*pi*fee.*fee;
C=exp(-1i*fee);
%C=exp(i*fee); %Time reversed
c=fsst(C); %clear R fee;
c=0.999*c/max(abs(c)); %normalize to abs-peak less 0.1%
C=fast(c);

L=C; %L(ko)=0; %zero Even freqs (0, 2fmin, 4fmin, ....)
R=C; %R(ke)=0;   %zero Odd freqs (fmin, 3fmin, 5fmin, ....)
l=fsst(L); r=fsst(R);

%Scale data to int32 2^31-1 range
G=.1;% set output
Gr=G; %Gr=0; %Overall gain
Gl=G; %Gl=0;
l=round(Gl*l*Nmax);
r=round(Gr*r*Nmax);
sprintf('Gains %d %d',l,r);

% l = zeros(1,length(l)); r = zeros(1,length(r)); %% send out nothing

%save the data in files stim0 and stim1

fd=fopen('data/stim0','w'); N=fwrite(fd,l,'int32'); %Scale to Nmax (32 bit integer)
fclose(fd);
fd=fopen('data/stim1','w'); N=fwrite(fd,r,'int32'); %Scale to Nmax (32 bit integer)
fclose(fd);
delete('./data/resp*') % delete any previously stored responses

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

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


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(['g7play -Q 0 -B ',num2str(NB),' -R0 data/stim0 -R1 data/stim1 -C0 data/resp0 -C1 data/resp1 -S %d &'],NFT); % chirp
    % MAC %command =sprintf(['g7play -Q 0 -B ',num2str(NB),' -R0 data/stim0 -R1 data/stim1 -C0 data/resp0 -C1 data/resp1 -S %d &'],NFT); % chirp
   
    %loopback 'none','server','i2s'
    %command =sprintf(['g7play -L none -Q 0 -B ',num2str(NB),' -F 1000 -M .1 -C0 data/resp0 -C1 data/resp1 -S %d &'],NFT); % sine wave
    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(1)
    %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(['g7play -Q 0 -B ',num2str(NB),' -R0 data/stim0 -R1 data/stim1 -U0 data/resp0 -U1 data/resp1 -S %d &'],NFT); % chirp
    
    % MAC %command =sprintf(['./g7play -Q 0 -B ',num2str(NB),' -R0 data/stim0 -R1 data/stim1 -U0 data/resp0 -U1 data/resp1 -S %d &'],NFT); % chirp
    
    %command =sprintf(['g7play -Q 0 -B ',num2str(NB),' -F 1000 -M .1 -U0 data/resp0 -U1 data/resp1 -S %d &'],NFT); % sine wave
    
    %loopback 'none','server','i2s'
    %command =sprintf(['g7play -L none -Q 0 -B ',num2str(NB),' -F 1000 -M .1 -U0 data/resp0 -U1 data/resp1 -S %d &'],NFT); % sine wave
    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(1)
        %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