make_ave_bf

function status=make_ave_bf(prestimsamps,prestimsamps2,lowpass,highpass,lowfiltord,highfiltord,infname,outfname)

Filter and baseline correct data in an average EGIS file.

This function accepts either 4, 6 or 8 input arguments. If the
last two arguments are omitted the function is run in 
interactive mode, and file names are obtained through the std
dialog boxes. 

The parameter prestimsamps1 & 2 are used to baseline correct through
a call to zeromeans (c.f.). The parameters lowpass and highpass
are used to construct the butterworth filter. If highpass is zero
than filtering occurs through a single lowpass filtering of order
20. Otherwise, data is bandpass filtered with a filter of order six,
and then refiltered with a lowpass filter of order 20 to improve
high-end rolloff.

Uses: ave_hdr_offsets_v, rd_egis_hdr_v, zeromeans, rd_onetr_allch

0001 function status=make_ave_bf(prestimsamps1,prestimsamps2,lowpass,highpass,lowfiltord,highfiltord,infname,outfname)
0002 %function status=make_ave_bf(prestimsamps,prestimsamps2,lowpass,highpass,lowfiltord,highfiltord,infname,outfname)
0003 %
0004 %Filter and baseline correct data in an average EGIS file.
0005 %
0006 %This function accepts either 4, 6 or 8 input arguments. If the
0007 %last two arguments are omitted the function is run in
0008 %interactive mode, and file names are obtained through the std
0009 %dialog boxes.
0010 %
0011 %The parameter prestimsamps1 & 2 are used to baseline correct through
0012 %a call to zeromeans (c.f.). The parameters lowpass and highpass
0013 %are used to construct the butterworth filter. If highpass is zero
0014 %than filtering occurs through a single lowpass filtering of order
0015 %20. Otherwise, data is bandpass filtered with a filter of order six,
0016 %and then refiltered with a lowpass filter of order 20 to improve
0017 %high-end rolloff.
0018 %
0019 %Uses: ave_hdr_offsets_v, rd_egis_hdr_v, zeromeans, rd_onetr_allch
0020 
0021 %Modification History
0022 %    9/95 Created by Ramesh Srinivasan
0023 %    10/95 Rewrote user interface and added high pas filtering - B. Rakitin
0024 %    10/31/95 Bug fixed in file opening branch - B. Rakitin
0025 %   11/17/95 reordered filters in bad-pass case - RS
0026 %    2/26/96 added pass params for filt orders - BCR
0027 
0028 status = -1;
0029 
0030 %Check number of input arguments
0031 if ~((nargin == 4)|(nargin == 6)|(nargin == 8))
0032     error('Number of input arguments must be either 4, 6, or 8.');
0033 end
0034 %Initialize fids
0035 srcfid = -1;
0036 destfid = -1;
0037 %First try batch mode
0038 if nargin == 8
0039     srcfid = fopen(infname, 'r');
0040     if srcfid == -1
0041         error(['Could not open input file' infname '.']);
0042     end
0043     destfid = fopen(outfname, 'wb');
0044     if destfid == -1
0045         temp =fclose(srcfid);
0046         error(['Could not open output file ' outfname '.']);
0047     end
0048     if highfiltord==[]
0049         highfiltord=3;
0050     end
0051     if lowfiltord==[]
0052         lowfiltord=20;
0053     end
0054 %Otherwise run interactive mode
0055 else
0056     while srcfid == -1
0057         [srcfid, infname]=get_fid('r','*ave*', 'Open Average File:');
0058     end
0059     while destfid == -1
0060         [destfid, outfname]=put_fid('wb','new.ave_bf','Save New Average File As:');
0061     end
0062     if nargin==6
0063         if highfiltord==[]
0064             highfiltord=3;
0065         end
0066         if lowfiltord==[]
0067             lowfiltord=20;
0068         end
0069     else
0070         highfiltord=3;
0071         lowfiltord=20;
0072     end
0073 end
0074 
0075 %Call EGIS hdr index script
0076 ave_hdr_offsets_v;
0077 %Read in the EGIS header
0078 [fhdr,chdr,ename,czeros,cgains,cnames,fcom,ftext, coff]=rd_egis_hdr_v(srcfid);
0079 %Calculate nyquist frequency and compare to lowpass parameter
0080 nyquist=zeros(1,fhdr(NCells));
0081 for c=1:fhdr(NCells)
0082     nyquist(c) = chdr(c, SampRate)/2;
0083 end
0084 if any(lowpass > nyquist)
0085     temp=fclose(srcfid);
0086     temp=fclose(destfid);
0087     error(['The lowpass freq. exceeds nyquist freq., ' int2str(nyquist) ', in some cells.']);
0088 end
0089 
0090 
0091 % Copy source header to destination file
0092 frewind(srcfid);
0093 [temp, countin] = fread(srcfid, fhdr(LHeader), 'char');
0094 countout = fwrite(destfid, temp, 'char');
0095 if countin ~= countout
0096   error('Error copying header information');
0097 end
0098 
0099 %Begin looping through cells
0100 for c=1:fhdr(NCells)
0101     disp(['Processing cell: ' int2str(c)]);
0102     %Construct the filters
0103     if highpass ~= 0
0104         %filter_order = 3;
0105         cutoff = [(highpass/nyquist(c)) (lowpass/nyquist(c))];
0106         b_bandpass=[]; a_bandpass=[];
0107         [b_bandpass, a_bandpass] = butter(highfiltord, cutoff);
0108             if c == 1
0109                 freqz(b_bandpass,a_bandpass);
0110                 figure
0111             end;
0112                 
0113     end;
0114 %    construct low-pass filter
0115         %filter_order = 20;
0116         cutoff = lowpass/nyquist(c);
0117         b_lowpass=[]; a_lowpass=[];
0118         [b_lowpass,a_lowpass] = butter(lowfiltord,cutoff);
0119             if c == 1
0120                 freqz(b_lowpass,a_lowpass);
0121             end;
0122 
0123     %Loop through subject averages
0124     for t=1:chdr(c,NObs)
0125         trialdata = rd_onetr_allch(srcfid, coff(c), t, fhdr(NChan), chdr(c, NPoints));
0126         for ch=1:fhdr(NChan)
0127                 filtdata(:,ch) = filtfilt(b_lowpass, a_lowpass, trialdata(:,ch));
0128         end;
0129     
0130         if highpass ~= 0
0131             for ch=1:fhdr(NChan)
0132                 filtdata(:,ch) = filtfilt(b_bandpass, a_bandpass, filtdata(:,ch));
0133             end;
0134         end;
0135         bfiltdata= zeros(size(filtdata,1),size(filtdata,2));
0136         bfiltdata = zeromean(filtdata,prestimsamps1,prestimsamps2);
0137         fwrite(destfid, bfiltdata', 'int16');
0138     end     %for t=1:chdr(c,NObs)
0139 end     %for c=1:fhdr(NCells)
0140 
0141 disp('Finished running make_ave_bf.');
0142 temp=fclose(srcfid);
0143 temp=fclose(destfid);
0144 status = 1;