%change lines 3,4,5,6,27,30,46-50 according to test
load('SI_sinusoidal_data.mat')
fulldata = Sinusoid64min_SI;  %define matrix to use
Period = '64 minute period';  %string of period length
Fs=1;   %sampling frequency (samples per minute)

P=65.28;    %pumping period length (found empirically)
Freq=1/P;   %pumping frequency
[N,m] = size(fulldata); %length,N and width, m of data
T=1/Fs;         %data collection time (minutes)
t=(0:N-1)*T;    % time axis in mins

Fstop = Freq*.5;       % Stopband Frequency
Fpass = Freq*.9;       % Passband Frequency
Astop =10;          % Stopband Attenuation (dB)
Apass = .1;           % Passband Ripple (dB)
match = 'stopband';  % Band to match exactly

% Construct an FDESIGN object and call its BUTTER method.
d = designfilt('highpassiir','StopbandFrequency',Fstop ,...
  'PassbandFrequency',Fpass,'StopbandAttenuation',Astop, ...
  'PassbandRipple',Apass,'SampleRate',Fs,'DesignMethod','butter');

clip = floor(Fs*2*P):floor(Fs*6*P); %clip first two Ps off and truncate at 8th
N_clip = length(clip); %length of clip
tc = (0:N_clip-1)*T; %time axis for clipped data
Nfft = N_clip-1; %FFT length for clipped data
Aph64min = zeros(4,m-2); %create blank matrix to store amp and phase
Bin_res = Fs/Nfft; %bin resolution
Hz=(0:Nfft-1)*(Fs/Nfft); %frequency axis
peak = 5; % bin number corresponding to peak at pumping freq
name = fulldata.Properties.VariableNames; %variable name for chart titles


for i = 2:m;
    filename = cell2mat(name(i));
    filtdata(:,i-1) = filtfilt(d,fulldata{:,i});  %filter the response
    data(:,i-1) = filtdata(clip,i-1);
    signal_mag = abs(fft(data(:,i-1),Nfft))/(.5*Nfft);
    signal_phase1 = unwrap(angle(fft(fulldata{:,i})));
    signal_phasewrap = unwrap(angle(fft(data(:,i-1),Nfft)));
    signal_phase = angle(fft(data(:,i-1),Nfft));
    phase(1,i-1) = signal_phase(peak);
    phase(2,i-1) = signal_phasewrap(peak);
    [M,L] = max(signal_mag);
    Aph64min(1,i-1) = L;
    Aph64min(2,i-1) = M;
    Aph64min(3,i-1) = signal_phase(peak);
    figure; plot(tc,data);xlabel('minutes');ylabel('head (m)');title([name(i),Period 'filtered time domain response'])
    figure; plot(Hz,signal_mag);xlabel('frequency (1/minute)');ylabel('amplitude (m)');title([name(i),Period 'FFT single sided amplitude spectrum']);...
     text(Freq,1.05*max(signal_mag),'\downarrow pumping frequency'); axis([0,Fs/2,0,1.1*max(signal_mag)])
    end

