diff --git a/Filter0.m b/Filter0.m new file mode 100644 index 0000000..16be29d --- /dev/null +++ b/Filter0.m @@ -0,0 +1,26 @@ +% y = Filter0(b, x) +% +% filters x with a fir filter so it has zero phase, i.e. shifts the +% filtered signal to the right half of the length of b. +% +% for now it zero pads the original signal +% later it might also do reflecton boundary conditions. +% +% be careful about the order of b! + + +function y = Filter0(b, x) + +if size(x,1) == 1 + x = x(:); +end + +if mod(length(b),2)~=1 + error('filter order should be odd'); +end + +shift = (length(b)-1)/2; + +[y0 z] = filter(b,1,x); + +y = [y0(shift+1:end,:) ; z(1:shift,:)]; diff --git a/SinoFilt.m b/SinoFilt.m new file mode 100644 index 0000000..4bf5ee6 --- /dev/null +++ b/SinoFilt.m @@ -0,0 +1,19 @@ +% SinoFilt, modified Matt Valley 110209 +% +% Usage: [Filtered,FiltAmp] = SinoFilt(eeg,sampFreq,filttype); +% filttype will be: [lowfreq highfreq] +% This will return the mean of the input +% channels filtered for the filttype + +function [Filtered,Filtamp] = SinoFilt(eeg,sampFreq,filttype,filtorder); + +Nyquist = sampFreq/2; +MyFilt=fir1(filtorder,filttype/Nyquist); + +%fprintf('\nFiltering...\n');% +Filtered = Filter0(MyFilt,eeg); + +if (nargout>1) + FiltHilb = hilbert(Filtered); + FiltAmp = abs(FiltHilb); +end diff --git a/spectral analysis/parse_xchan.asv b/spectral analysis/parse_xchan.asv new file mode 100644 index 0000000..43fa319 --- /dev/null +++ b/spectral analysis/parse_xchan.asv @@ -0,0 +1,71 @@ +clear all; + +% Define Global VARS +eventcodes = [1 2 3]; +srate = 3051.76; %Hz +winsize = .5; +brthindx = [-10:1:20];; +maxfreq = 120; +dead_chans = [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]; + %1 2 3 4 5 6 7 8 9 1 1 2 3 4 5 6 7 8 9 1 1 2 3 4 5 6 7 8 9 1 1 2 +%sigtype = 0; %0 gives mean power, 1 gives max power + +base_mode = 1; +% 0 = normalizes to average baseline power in each band individaully +% 1 = normalizes to average baseline power in each band individaully, and sets all baseline specs to zero +% 2 = normalizes to max baseline power of all gamma freq bands +% 3 = normalizes to max baseline power averaged over gamma freq bands +% 4 = normalizes to average baseline power in each band individaully, and also adjusts for band dynamic range +% 5 = normalizes to average baseline power in each band individaully, and also adjusts for band dynamic range after db normalization + +% upload the breath file - "breaths" in workspace +[datafile, pathname] = uigetfile(... + '*.mat',... + 'Please pick breath file'); +for n = 1:length(datafile); + cd(pathname); + load(datafile); +end +% upload the event file - "events" in workspace +[datafile, pathname] = uigetfile(... + '*.mat',... + 'Please pick event file'); +for n = 1:length(datafile); + cd(pathname); + load(datafile); +end + +% sequentially open all channels +[datafile, pathname] = = uipickfiles('\.mat$', '*.mat', 'MAT-files' }) +cd(pathname); +datafile = datafile([2:32 1]); %fix uigetfile bug, 32chan + +for odor = eventcodes; + for n = 1:length(datafile); %n to number of channels + if dead_chans(n)==1 + continue + end + load(datafile{n}); % name of array must be "wave" + for x = 1:length(brthindx) %x to num of breaths + wave_segs(:,:,x,n) = parsechans(wave,events,breaths,srate,odor,brthindx(x),winsize); + [S, t, f] = pmtm_cust(squeeze(wave_segs(:,:,x,n)),srate,maxfreq); + spec(:,:,:,x) = S; clear S; %spec = (time, freq, trial, breath) + disp('breath'); disp(x); + end + g_freqs = find(f>50 & f<100); + disp('site'); disp(n); + spec_all(:,:,:,:,n) = spec; + [spec_norm(:,:,:,:,n) aveallgamma(:,:,n)] = pmtmprocess_111709(spec,f,brthindx,base_mode,g_freqs); + [sig_breaths(:,:,n),sig_vals(:,:,n),cis(:,:,:,n),all_spec(:,:,:,n)] = test_breathsig(spec_norm(:,:,:,:,n),brthindx,g_freqs); + end + aveallgamma_allodors(:,:,:,odor) = aveallgamma; + sig_breaths_allodors(:,:,:,odor) = sig_breaths; + save(['spec_odor' num2str(odor)], 'spec_all'); + clear spec_all; + save(['spec_norm_odor' num2str(odor)], 'spec_norm'); + clear spec_norm; + save(['wave_segs' num2str(odor)], 'wave_segs'); + save('f','f'); +end +save(['aveallgamma_allodors'], 'aveallgamma_allodors'); +save(['sig_breaths_allodors'], 'sig_breaths_allodors'); \ No newline at end of file diff --git a/spectral analysis/parse_xchan.m b/spectral analysis/parse_xchan.m index e4990df..4ff260a 100644 --- a/spectral analysis/parse_xchan.m +++ b/spectral analysis/parse_xchan.m @@ -4,8 +4,9 @@ eventcodes = [1 2 3]; srate = 3051.76; %Hz winsize = .5; -brthindx = [-10:1:20];; +brthindx = [-10:1:20]; maxfreq = 120; +gamma = [50 100]; %define gamma frequency window (Hz) dead_chans = [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]; %1 2 3 4 5 6 7 8 9 1 1 2 3 4 5 6 7 8 9 1 1 2 3 4 5 6 7 8 9 1 1 2 %sigtype = 0; %0 gives mean power, 1 gives max power @@ -35,7 +36,7 @@ load(datafile); end -% sequentially open all channels +% sequentially import all channels [pathname] = uipickfiles('refilter', '\.mat$', 'type', {'*.mat', 'MAT-files'},... 'prompt', 'Select all .mat files from MEA channels', 'output', 'cell'); @@ -45,14 +46,14 @@ if dead_chans(n)==1 continue end - load(pathname{n}); % name of array must be "wave" + load(pathname{n}); % name of imported data must be "wave" for x = 1:length(brthindx) %x to num of breaths wave_segs(:,:,x,n) = parsechans(wave,events,breaths,srate,odor,brthindx(x),winsize,eventcodes); [S, t, f] = pmtm_cust(squeeze(wave_segs(:,:,x,n)),srate,maxfreq); spec(:,:,:,x) = S; clear S; %spec = (time, freq, trial, breath) %disp('breath'); disp(x); end - g_freqs = find(f>50 & f<100); + g_freqs = find(f>gamma(1) & f