function [s] = miso_filter_int(s,y)
%miso_filter_int inputs a filter and a structure array of data sets y, applies filter to data, and
%outputs a structure with fields: ppos(signal frequ spectrum),perr(cost
%function frequ spectrum),pest(signal estimate frequency
%spectrum),f(frequency),target(signal),est_darm(noise estimate),t(time).
%data file for which filter has been calculated is s (obtained using miso_filter). 
%y consists of data structures which will be filtered using
%filter from s. Then the power spectrum of the difference between signal and filtered-data is
%graphed for all the data files of y for comparison too see how well filter performs
%over time. Note if you want to create a y, take z1,z2,z3,etc. structures
%obtained by using get_data, and let y=[z1;z2;z3]. Note that s also contains
%information regarding how well the filter performed, like a noise
%estimate. This info is included in the graph at the end of miso_filter_int. 


%reading in s
wiener_input_data=s.noise;
h=s.filter;
srate=s.sample_rate;
target=s.signal;
rat=s.rat;
N=s.order;
duration=s.duration;
t0=s.start;
t=s.time;
est_darm=s.est_noise;
ppos=s.f_signal;
pest=s.f_est;
perr=s.f_err;
f=s.freq;

%Saving memory
clear s

%determining how many data stuctures there are to be filered.
x = size(y);
x = x(1);

% Loop over data structures
for ii= 1:x
 z = y(ii,:);

%loop over channels and downsampling and detrending data to be filtered and
%noise data
 disp('Downsampling and detrending')
 for n = 1:length(z)
   z(n).data = decimate(z(n).data(:), z(n).rate/srate);
   z(n).data = detrend(z(n).data);
 end
 
% CAT the input data columns into a matrix
 wiener_input_data_temp = [];
 for gg = 1:length(z)-1
   wiener_input_data_temp = [wiener_input_data_temp z(gg).data];
 end
 target_temp = z(end).data;

% DO Time domain FIR filtering using the Wiener coeffs for each time interval
 disp('Applying filter')
 est_darm_temp = 0;
 [a, wl] = size(z);

 % Loop over channels to filter signal 
 for n = 1:wl-1
   est_darm_temp = est_darm_temp +...
             filter(h(1+N*(n-1):n*N), 1, wiener_input_data_temp(:,n));
 end

 %making time axis to plot target and est_darm with later
 t_temp= z(end).start + ((0:length(target_temp)-1)/srate)';
  
 %creating vector of signals and estimates of signals for each time interval
 est_darm=[est_darm,est_darm_temp];
 target=[target,target_temp];
 t=[t,t_temp];
 
 % PWELCH to prove that something good is happening
  nfft = srate*16;
  [ppos_temp,f_temp] = pwelch(target_temp,hanning(nfft),nfft/2,nfft,srate);
  [perr_temp,f_temp] = pwelch(target_temp - est_darm_temp,hanning(nfft),nfft/2,nfft,srate);
  [pest_temp,f_temp] = pwelch(est_darm_temp,hanning(nfft),nfft/2,nfft,srate);

  ppos_temp = sqrt(ppos_temp);
  perr_temp = sqrt(perr_temp);
  pest_temp = sqrt(pest_temp);
  
 %making vectors of ppos,perr,pest for all time intervals 
  ppos=[ppos,ppos_temp];
  pest=[pest,pest_temp];
  perr=[perr,perr_temp];
  f=[f,f_temp];

% clearing info because of loop
  clear target_temp est_darm_temp ppos_temp perr_temp pest_temp z ...
      wiener_input_data_temp wl nfft f_temp
 
end

%saving memory
clear y

s=struct('ppos',ppos,'pest',pest,'perr',perr,'f',f,'target',target,'est_darm',est_darm,'t',t)

 % PLOT difference in first 5 intervals - - - - - - - - - - - - - - - 
if x>4 
 loglog(f(:,1),rms(f(:,1),perr(:,1)),'r',...
        f(:,2),rms(f(:,2),perr(:,2)),'m',...
        f(:,3),rms(f(:,3),perr(:,3)),'b',...
        f(:,4),rms(f(:,4),perr(:,4)),'g',...
        f(:,5),rms(f(:,5),perr(:,5)),'c',...
        f(:,6),rms(f(:,5),ppos(:,5)),'k')    

 grid
 xlabel('Frequency [Hz]')
 ylabel('cts/rHz')
 axis tight
 title('MISO Wiener Filter Time')
 legend('residual 0',...
       'residual 1',...    
       'residual 2',... 
       'residual 3',...      
       'residual 4',... 
       'target 4',...
       'Location','SouthWest') 
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
end


%plot(t,target,'red',t,est_darm,'green')