function [s] = miso_filter_lev(N,srate,rat,z)
%MISO_FILTER_LEV(N,srate,z) uses miso_firlev to get levinson
%   FIR Wiener filter of order N-1, using impulse response of 
%   N/srate. z is a structure gotten from the get_data function. 
%   z(end) is the signal which is filtered using z(i) for all i.
%   'rat' is the fraction of z which will be put into filter
%   funtion. The data from z is downsampled using srate and 
%   detrended. Let rat=1. I don't have that part working yet.


rat=1;

disp('Downsampling and detrending')
for n = 1:length(z)
  z(n).data=z(n).data';
  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 = [];
for gg = 1:length(z)-1
  wiener_input_data = [wiener_input_data z(gg).data];
end
target = z(end).data;


% Getting start time and duration and name
t0 = z.start;
duration = z.duration;
name=z.name


% SAVE Memory
clear z 


%Calc filter
disp('Calculating Wiener filter')
tic
h = miso_firlev(N-1,wiener_input_data,target);
toc

% DO Time domain FIR filtering using the Wiener coeffs
disp('Applying filter')
est_darm = 0;
[a, wl] = size(wiener_input_data);
for n = 1:wl
   est_darm = est_darm +...
             filter(h(1+N*(n-1):n*N), 1, wiener_input_data(:,n));
end

% PWELCH to prove that something good is happening
nfft = srate * 16;
[ppos,f] = pwelch(target,hanning(nfft),nfft/2,nfft,srate);
[perr,f] = pwelch(target - est_darm,hanning(nfft),nfft/2,nfft,srate);
[pest,f] = pwelch(est_darm,hanning(nfft),nfft/2,nfft,srate);
ppos = sqrt(ppos);
perr = sqrt(perr);
pest = sqrt(pest);

% PLOT stuff just to show off - - - - - - - - - - - - - - - 
semilogy(f,ppos,'r',...
       f,rms(f,ppos),'r--',...
       f,perr,'k',...
       f,rms(f,perr),'k--',...
       f,pest,'b',...
       f,rms(f,pest),'b--')

grid
xlabel('Frequency [Hz]')
ylabel('cts/rHz')
axis tight
%axis([0.1 60 0.01 200])
title('MISO Wiener Filter')
nam=name(end,:)
legend([nam],...
       [nam ' (rms)'],...
       'residual',...
       'residual (rms)',...
       'est',...
       'est(rms)',...
       'Location','SouthWest')
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 

% making a time axis to plot data with later
t=t0 + ((0:length(target)-1)/srate)';
%plot(t,target,'red',t,est_darm,'blue',t,target-est_darm,'green')
%plot(t,target,'red',t,target-est_darm,'green',t,est_darm,'m')


% writing info about structure s below
info=['start(start time),duration,time(time axis), filter(output from'] ;
info=[info ' miso_firwiener), noise(inputs of noise for miso),signal(original signal'];
info=[info ' which was being filtered),est_signal(signal after being filtered by'];
info=[info ' using miso),f_signal(power spectrum of signal),f_est(power spectrum'];
info=[info ' of est_signal), f_err(power spectrum of signal-est_signal), order(order'];
info=[info ' of filter), sample rate, rat(fraction of signal used in calculating'];
info=[info ' filter)'];

% creating structure to save
s=struct('name',name,'start',t0,'duration',duration,'time',t,...
         'filter',h,'noise',wiener_input_data,'signal',target,...
         'est_noise',est_darm,'f_signal',ppos,...
         'f_est',pest,'f_err',perr,'freq',f,'order',N,...
         'sample_rate',srate,'rat',rat,'info',info);

% Making a name which includes start time and duration
name = ['w_'  nam '_' num2str(t0) '_' num2str(duration) '_' num2str(srate) '_' ...
        num2str(N)]


save temp_w s

disp('data saved as temp_w. please save as name(above)')