xcorrelation work - continuing to try to make it functional

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,:)];

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

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

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<gamma(2));
         disp('site'); disp(n);
         spec_all(:,:,:,:,n) = spec;
         [spec_norm(:,:,:,:,n) aveallgamma(:,:,n)] = pmtmprocess(spec,f,brthindx,base_mode,g_freqs);

spectral analysis/parsechans.m

@@ -1,4 +1,4 @@
-function [wave_segs] = parsechans(rawwave,events,breaths,srate,odor,brth_num,winsize,eventcodes)
+function [wave_segs] = parsechans(wave,events,breaths,srate,odor,brth_num,winsize,eventcodes)
 
 %INPUTS:
 % -- "events" matrix containing vectors with event time(in sec) 1st dim,
@@ -24,7 +24,7 @@
     sel_events = sel_events(1:(length(sel_events)-1)); %cut out last event because the experiment is often ended before a sufficient number of breaths are recorded after this event
 end
 
-for a = 1:length(eventcodes);
+for a = eventcodes;
     b = find(tdt_allevents(:,2) == a);
     c(a) = length(b);
     clear b;
@@ -50,6 +50,6 @@
     if i==1;
         win_final = winend-winstart; %a bug occurs rarely where the int32 of the window creates a size mismatch between trials.  "winfinal" is the template size to maintain over all trials.
     end  
-    wave_segs(:,i) = rawwave(winstart:(winstart+win_final));
+    wave_segs(:,i) = wave(winstart:(winstart+win_final));
 end
 end

xcorrelation/CorrVDist.m

@@ -1,4 +1,4 @@
-function[corrs_distances]= CorrVDist(Corr_Per_Breath,coordinates,FirstChannel,LastChannel,ReferenceChannel,DeadChannel,Breath)
+function[corrs_distances]= CorrVDist(Corr_Per_Breath,coordinates,FirstChannel,LastChannel,ReferenceChannel,DeadChannel)
 
 %this is the function that plots correlation vs distance per specific
 %reference channel. Final version. It uses the corr _Per_ Breath output
@@ -26,7 +26,7 @@
 
 [distance_Chagit]=distcalc_perchan(coordinates,FirstChannel, LastChannel, ReferenceChannel,DeadChannel);    %calculates distance with respect to reference channel
 unique_final_distance=unique(distance_Chagit);          %unique all the distances
-unique_corrs = squeeze(Corr_Per_Breath(ReferenceChannel,:,Breath));
+unique_corrs = squeeze(Corr_Per_Breath(ReferenceChannel,:));
 
 for i=1:length(unique_final_distance)
     a=find(unique_final_distance(i)==distance_Chagit);                            %find where distance equal all the successive unique distances

xcorrelation/scatter_corrs.m

@@ -2,39 +2,46 @@
 %2    32    30    32
 events = size(Corr_Per_Breath,3);
 breaths = size(Corr_Per_Breath,4);
+basebreaths = 1:11;
+odorbreaths = 12:14;
+postbreaths = 15:30;
 % Corr_Per_Breath = [channels, channels, events, breaths]
 
-for event = 12:17;
+for e = 1:events;
     for chan=1:32;
-        for breath = 1:breaths;
-            [corrs_distances]= CorrVDist(squeeze(Corr_Per_Breath(:,:,events,:)),coordinates,1,32,chan,33,breath);
-            allcorrs(:,1:size(corrs_distances,2),breath,chan) = corrs_distances;
+        for b = 1:breaths;
+            [corrs_distances]= CorrVDist(Corr_Per_Breath(:,:,e,b),coordinates,1,32,chan,33);
+            allcorrs(:,1:size(corrs_distances,2),b,chan) = corrs_distances;
         end
     end
 end
 
-for breath = 1:10;
-    for chan = 1:32;
-        %i = (allcorrs(1,:,breath,chan),allcorrs(2,:,breath,chan),8);
-            %dist_vals = find(
-        %scatter(allcorrs(1,:,breath,chan),allcorrs(2,:,breath,chan),8);
-        distances = allcorrs(2,:,breath,chan);
-        plot(allcorrs(1,:,breath,chan), distances);
-        hold on;
-    end
-end
-
-mean_corrs = squeeze(mean(mean(allcorrs(1,:,:,:),3),4));
+% for breath = 1:10;
+%     for chan = 1:32;
+%         %i = (allcorrs(1,:,breath,chan),allcorrs(2,:,breath,chan),8);
+%             %dist_vals = find(
+%         %scatter(allcorrs(1,:,breath,chan),allcorrs(2,:,breath,chan),8);
+%         distances = allcorrs(2,:,breath,chan);
+%         plot(allcorrs(1,:,breath,chan), distances);
+%         hold on;
+%     end
+% end
+% 
+mean_corrs_base = squeeze(mean(mean(allcorrs(1,:,basebreaths,:),3),4));
+mean_corrs_odor = squeeze(mean(mean(allcorrs(1,:,odorbreaths,:),3),4));
+mean_corrs_post = squeeze(mean(mean(allcorrs(1,:,postbreaths,:),3),4));
 distances = allcorrs(2,:,1,1);
 
-for i = 1:size(allcorrs,2)
-    a = squeeze(allcorrs(1,i,:,:));
-    b(:,i) = a(:);
-end
-stdev_corrs = std(b,1);%./sqrt(size(b,1));
+% for i = 1:size(allcorrs,2)
+%     a = squeeze(allcorrs(1,i,:,:));
+%     b(:,i) = a(:);
+% end
+% stdev_corrs = std(b,1);%./sqrt(size(b,1));
 
-plot(distances,mean_corrs);
+plot(distances,mean_corrs_base, 'k');
 hold on;
-plot(distances,mean_corrs+stdev_corrs,'r');
+plot(distances,mean_corrs_odor, 'r');
+%plot(distances,mean_corrs+stdev_corrs,'r');
 hold on;
-plot(distances,mean_corrs-stdev_corrs,'r');
+plot(distances,mean_corrs_post, 'g');
+%plot(distances,mean_corrs-stdev_corrs,'r');