%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Creates a bode plot for a circuit driven
%%% by the audio cable. Input is hooked to ai0
%%% and output to ai0.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function bodeplot_exp()
    h1 = figure(1); clf; h2  = figure(2); clf;
    set(h1,'Position',[100   300   560   420]);
    set(h2,'Position',[700  300   560   420]);


    freq = logspace(0,3,20);   %% set up frequency sweep in logspace


    %%% parameters for the experiment
    sample_rate = 24000;

    if (sample_rate/max(freq) < 4)
        fprintf('Cannot give reliable results at freq of %f\n at sample rate of %f\n',max(freq),sample_rate)
        error('reset freq range and try again')
    end


    %%% set up the DAQ session
    d = daq.getDevices;  d =d(1);  %%% Get the devices you have plugged in, assume only one device
    s = daq.createSession('ni');                       %%% start the session
    s.addAnalogInputChannel(d.ID, 'ai0', 'Voltage');   %%% Add channel ai0
    s.Channels(1).Range = [-5 5];                      %%% set input range for ai0
    s.addAnalogInputChannel(d.ID, 'ai1', 'Voltage');   %%% Add channel ai1
    s.Channels(2).Range = [-5 5];                      %%% set input range for ai1
    s.Rate = sample_rate;           %% set the sample rate



    for i = 1:length(freq)
        duration = max(2,ceil(10/freq(i)));   %% how long to collect data
        s.DurationInSeconds = duration;
        lh = s.addlistener('DataAvailable',@(src, event)processData(src, event, freq(i)));
        s.NotifyWhenDataAvailableExceeds = s.NumberOfScans;

        s.startBackground;

        %%% set up signal to go out on audio -
        Fs = 44100;
        t = linspace(0,duration,Fs*duration )';
        y1 = sin(2*pi*freq(i)*t);

        soundsc([y1 y1 ],Fs);  %% send the audio signal

        s.wait();

        delete(lh);


    end

end


function processData(src, event,freq)

    N = length(event.TimeStamps);

    %%% extract the last half of signal to ignore transient
    input  = event.Data(floor(N/2):end,1);
    output = event.Data(floor(N/2):end,2);
    t = event.TimeStamps(floor(N/2):end);
    t = t-t(1);

    figure(1);
    subplot(2,1,1)
    plot(t,input); hold on; plot(t,output,'r'); hold off
    legend('input' ,'output')
    xlabel('time (s)'); ylabel('Voltage');
    axis([0 1*5/freq -max(abs(input))*1.1 max(abs(input))*1.1])
    subplot(2,1,2)
    plot(t,output,'r'); hold off
    legend('output')
    xlabel('time (s)'); ylabel('Voltage');
    axis([0 1*5/freq  -max(abs(output))*1.1 max(abs(output))*1.1])
    drawnow;


    %%% compute the sine and cosine components at the driving freq
    A_in =  2*trapz(t,input.*sin(2*pi*freq*t))/(t(end)-t(1));
    B_in =  2*trapz(t,input.*cos(2*pi*freq*t))/(t(end)-t(1));
    A_out = 2*trapz(t,output.*sin(2*pi*freq*t))/(t(end)-t(1));
    B_out = 2*trapz(t,output.*cos(2*pi*freq*t))/(t(end)-t(1));

    %%% Compute the amplitude and phase
    Amp_out = sqrt(A_out^2+B_out^2);
    Amp_in = sqrt(A_in^2+B_in^2);
    Amp = Amp_out/Amp_in;
    Phase = angle(complex(A_in,B_in)) - angle(complex(A_out,B_out));
    Phase = Phase*180/pi;
    if (Phase<0)
        Phase = Phase+360;
    end

    %%% create the Bode plot
    figure(2);
    subplot(2,1,1)
    loglog(freq,Amp,'.'); hold on;
    xlabel('Frequency (Hz)')
    ylabel('Amplitude out/Amplitude in')
    subplot(2,1,2)
    semilogx(freq,Phase,'.'); hold on;
    xlabel('Frequency (Hz)')
    ylabel('Phase (degrees)')
    drawnow

    fprintf('Freq %10.2f Amp in %9.6f  Amp out %9.6f Phase %7.2f \n',...
        freq,Amp_in,Amp_out,Phase);

end