%last modified 5/31/06
close all;
clear all;
clc;

%define string arrays for easy file accessability
names = ['NEWWPC_MODALTEST01.txt';'NEWWPC_MODALTEST02.txt'];%'hit31.txt'];%;'test4.txt';'test5.txt';'test6.txt'];
[k,m] = size(names);     %obtain number of rows 'k'

%initialize transfer function values
Sum_tf = 0;

%calculate transfer function at the tip
for i = 1:1:k
    data = load(names(i,:));
    t   = data(:,1);
    hv   = data(:,3);
    av   = data(:,2);



    acal    = 1/0.979e-3;   %acceleration calibration constant, m/s^2/V
    hcal    = 2.330005368921638e+002;   %impact calibration constant, N/V
    a   = av*acal;          %acceleration, m/s^2
    h   = hv*hcal;          %impact, N
    
    
    figure(11);
    subplot(2,1,1), plot(t,h) 
    subplot(2,1,2), plot(t,a)
    subplot(2,1,1); zoom on;
    pause;
    [tval,dummy] = ginput(2);      
    h = 1./(1+exp(-600*(t-tval(1))))*1./(1+exp(600*(t-tval(2)))).*h; 
    hold on; 
    subplot(2,1,1), plot(t,h,'r')
    subplot(2,1,2); zoom on;
    pause(1);
    [tval,dummy] = ginput(2);
    a = 1./(1+exp(-600*(t-tval(1))))*1./(1+exp(600*(t-tval(2)))).*a;
    hold on; 
    subplot(2,1,2), plot(t,a,'r')
    pause(2)
    close(11)
    
    

    % TAKE THE FOURIER TRANSFORM OF THE SIGNALS
    npts    = length(h);
    hf      = 2*fft(h)/npts;    % this is the raw fourier transfer
    hf(1)   = hf(1)/2;          % correction for the zero frequency
    fs      = 1/(t(2)-t(1));    % sample rate in hertz
    fhz     = fs*(t-t(1))/(max(t)-t(1));

    af      = 2*fft(a)/npts;   % raw fourier transform for acceleration
    af(1)   = af(1)/2;          % correction for the zero frequency

    w       = 2*pi*fhz;
    w2      = w.*w;

    xf      = -af./w2;  %displacement of beam, m
    tf      = xf./hf;   %displacement to force transfer function, m/N
    
    Sum_tf = Sum_tf+tf; %sum tip transfer functions for each test
end
Tf = Sum_tf/k;  %Average tip transfer function
figure(3); subplot(2,1,1), plot(fhz,real(Tf));
xlim([100,7000])
title('Experimental Response')
xlabel('Freq (Hz)')
ylabel('Real (m/N)')
subplot(2,1,2), plot(fhz,imag(Tf));
xlim([100,7000])
xlabel('Freq (Hz)')
ylabel('Imag (m/N)')

HH(1,:) = real(Tf);
HH(2,:) = imag(Tf);

save file3 Tf fhz HH