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)')