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