reduce3

% Data reduction routine for MMT elevation open-loop data
% generally, we have the following:
%   Sample rate = 1000Hz
%   xPC Target version 2.7.1 (i.e. Matlab v7.04.1)
% The open-loop code logs the input outdata as:
%   1 => Input voltage
%   2 => Input frequency
%   3 => Absolute encoder, degrees
%   4 => Drive arc encoder, counts
%   5 => Drive arc encoder, counts
% Inputs are volts to the current-source linear amplifier, counts are
% interpolated 50X, and degrees calculated from the meshed raw absolute
% encoder outputs. Interpolation @50X gives 364007 counts/degree.
%
%       MAT-file for this script: El_open_compare121206.mat
%
%       additional data path == \\Mmto\dclark\MMT Servo Data\Comparative
%       file: MMT_el_open33005.mat
%

%generalize the i/o variables and tfestimate routine
    input9 = dataf9(:,1);
    output1 = dataf9(:,4) * 1/364007;
    output2 = dataf9(:,5) * 1/364007;

    freqs9 = dataf9(:,2);

    input5 = dataf5(:,1);
    output3 = dataf5(:,4) * 1/364007;
    output4 = dataf5(:,5) * 1/364007;

    freqs5 = dataf5(:,2);

%imported older data file for comparative purposes, f/5 hecto was on.
    ostart = 1.5e5;
    oldfreqs = olddata(:,2);
    outold = olddata(ostart:end,5) * 1/364007; %W tape @ 50X, clip "slide" at start
    inold = olddata(ostart:end,1);


%subtract any DC offset using the last 10K samples to generate the average
len = length(output1) - 1e4;
mean_1 = mean(output1(len:end));
output1 = output1 - mean_1;

len = length(output2) - 1e4;
mean_2 = mean(output2(len:end));
output2= output2 - mean_2;

len = length(output3) - 1e4;
mean_3 = mean(output3(len:end));
output3 = output3 - mean_3;

len = length(output4) - 1e4;
mean_4 = mean(output4(len:end));
output4 = output4 - mean_4;

len = length(outold) - 1e4;
mean_5 = mean(outold(len:end));
outold = outold - mean_5;

%tfe setup variable
win = 1024; fftlen = 16384; rate = 1000;

%calculate the current new transfer function estimate
    [Tf9a,Ff9a] =  tfestimate(input9,output1,win,[],fftlen,rate);
    Mag_f9a = 20 * log10(abs(Tf9a));
    Phs_f9a = unwrap(angle(Tf9a)) * 180/pi;

    [Tf9b,Ff9b] =  tfestimate(input9,output2,win,[],fftlen,rate);
    Mag_f9b = 20 * log10(abs(Tf9b));
    Phs_f9b = unwrap(angle(Tf9b)) * 180/pi;

    [Tf5a,Ff5a] =  tfestimate(input5,output3,win,[],fftlen,rate);
    Mag_f5a = 20 * log10(abs(Tf5a));
    Phs_f5a = unwrap(angle(Tf5a)) * 180/pi;

    [Tf5b,Ff5b] =  tfestimate(input5,output4,win,[],fftlen,rate);
    Mag_f5b = 20 * log10(abs(Tf5b));
    Phs_f5b = unwrap(angle(Tf5b)) * 180/pi;

    [Tf5old,Ffold] =  tfestimate(inold,outold,win,[],fftlen,rate);
    Mag_f5old = 20 * log10(abs(Tf5old));
    Phs_f5old = unwrap(angle(Tf5old)) * 180/pi;

%whack things around a bit for ez display later, i.e.
%bring the frequency vector down below Fs/2 and only display the results
%within the frequencies of the original input signal
Finlo9 = min(freqs9(1:1000));
Finhi9 = max(freqs9);
Flo9 = min(find(Ff9 > Finlo));
Fhi9 = min(find(Ff9 > Finhi));

Finlo5 = min(freqs5(1:1000));
Finhi5 = max(freqs5);
Flo5 = min(find(Ff5 > Finlo));
Fhi5 = min(find(Ff5 > Finhi));

Finloold = min(oldfreqs(1:1000));
Finhiold = max(oldfreqs);
Floold = min(find(Ffold > Finloold));
Fhiold = min(find(Ffold > Finhiold));

%plot results

figure(1);

title 'El Axis Open-loop Response';
    semilogx(Ff5a(Flo5:Fhi5),Mag_f5a(Flo5:Fhi5),'b');
    hold on;
    semilogx(Ff9a(Flo9:Fhi9),Mag_f9a(Flo9:Fhi9),'r');
    semilogx(Ff5b(Flo5:Fhi5),Mag_f5b(Flo5:Fhi5),'-.b');
    semilogx(Ff9b(Flo9:Fhi9),Mag_f9b(Flo9:Fhi9),'-.r');
    semilogx(Ffold(Floold:Fhiold),Mag_f5old(Floold:Fhiold),'g');
ylabel 'Db mag, counts';
grid on;
legend 'f/5' 'f/9' 'f/5' 'f/9' '2005 hecto';

figure(2);

    semilogx(Ff5a(Flo5:Fhi5),Phs_f5a(Flo5:Fhi5),'b');
    hold on;
    semilogx(Ff9a(Flo9:Fhi9),(180 - Phs_f9a(Flo9:Fhi9)),'r');
    semilogx(Ff5b(Flo5:Fhi5),Phs_f5b(Flo5:Fhi5),'-.b');
    semilogx(Ff9b(Flo9:Fhi9),(180 - Phs_f9b(Flo9:Fhi9)),'-.r');
    semilogx(Ffold(Floold:Fhiold),Phs_f5b(Floold:Fhiold),'g');
ylabel 'Phase, degrees';
grid on;
legend 'f/5' 'f/9' 'f/5' 'f/9' '2005 hecto';