-
Notifications
You must be signed in to change notification settings - Fork 3
/
analyzingTravellingWaves.m
115 lines (93 loc) · 3.29 KB
/
analyzingTravellingWaves.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
function analyzingTravellingWaves(avec)
%avec comes from echoPred
if length(size(avec))==3
randTrial=ceil(size(avec,2)*rand(1));
avec=squeeze(avec(:,randTrial,:));
end
%few parameters to set
timeWindow=200;
bins=10;
x=1:timeWindow:(timeWindow*bins)+1;
a=1; %width pixel to compute the frequency
arrangedAvec=avec';
avecFFT=abs(fftshift(fft2(arrangedAvec)));
[m,n]=size(arrangedAvec); %numberPixel x and y
fx=(1/a)*((1:n)-mean(1:n));
fy=(1/a)*((1:m)-mean(1:m));
[thX1,thX2]=findingXboundaries(avecFFT,fx);
figure
subplot(2,1,1)
surf(arrangedAvec)
shading interp; view(0,90); axis tight; colorbar;
xlabel('time [ms]')
ylabel('layers')
title('Looking for travelling waves')
subplot(2,2,3)
surf(fx,fy,avecFFT)
title('2D FFT - zoomed xAxis')
xlabel('timeFreq')
ylabel('layersFreq')
shading interp; view(0,90); axis tight; colorbar;
set(gca,'xlim',[thX1 thX2],'clim',[min(min(avecFFT)) max(max(avecFFT))],'yscale','lin');
subplot(2,2,4)
avecFFT(avecFFT==0)=NaN;
upperPart=avecFFT(fy>0,(fx>0 & fx<thX2));
lowerPart=avecFFT(fy<0,(fx>0 & fx<thX2));
bar([1 2],[mean(upperPart(:)) mean(lowerPart(:))])
title('Quantifying the waves')
xlabel('upper right quadrant x>0, y>0 - lower right quadrant x>0, y<0')
ylabel('sum of the quadrant')
m=8;
n=timeWindow;
a=1;
fx=(1/a)*((1:n)-mean(1:n));
fy=(1/a)*((1:m)-mean(1:m));
yy=0.001;
figure
for ii=1:length(x)-1
tempAvec=arrangedAvec(:,x(ii):x(ii+1)-1);
tempAvecFFT=abs(fftshift(fft2(tempAvec)));
[thX1temp,thX2temp]=findingXboundaries(tempAvecFFT,fx);
subplot(3,length(x),ii)
surf(x(ii):x(ii+1)-1,1:8,tempAvec)
shading interp; view(0,90); axis tight;
set(gca,'clim',[min(min(tempAvecFFT)) max(max(tempAvec))],'yscale','lin');
tempAvecVect(:,:,ii)=tempAvec;
subplot(3,length(x),ii+length(x))
surf(fx,fy,tempAvecFFT)
shading interp; view(0,90); axis tight;
set(gca,'xlim',[thX1temp thX2temp],'clim',[min(min(tempAvecFFT)) max(max(tempAvecFFT))],'yscale','lin');
subplot(3,length(x),ii+length(x)*2)
tempAvecFFT(tempAvecFFT==0)=NaN;
upperPart=tempAvecFFT(fy>0,(fx>0 & fx<thX2temp));
lowerPart=tempAvecFFT(fy<0,(fx>0 & fx<thX2temp));
bar([1 2],[mean(upperPart(:)) mean(lowerPart(:))])
timeVect(ii,:)=[mean(upperPart(:)) mean(lowerPart(:))];
% set(gca,'ylim',[-0.05 0.05],'yscale','lin');
end
figure
for ii=1:length(x)-1
subplot(2,ceil(length(x)/2),ii)
surf(x(ii):x(ii+1)-1,1:8,tempAvecVect(:,:,ii))
shading interp; view(0,90); axis tight;
set(gca,'clim',[min(min(min(tempAvecVect))) max(max(max(tempAvecVect)))],'yscale','lin');
title(['2D maps - timeBin: ' int2str(ii)])
xlabel('time [ms]')
ylabel('layers')
end
subplot(2,ceil(length(x)/2),length(x))
plot(timeVect)
title('2D FFT - quantifying')
legend('upper quadrant','lower quadrant')
xlabel('time bins')
ylabel('sum for each quadrant')
end
function [thX1,thX2]=findingXboundaries(avecFFT,fx)
%this is just a makiavelic way to find the boundaries of the interesting frequencies (i.e. thX1,thX2) ##to improve
sumFFT=sum(avecFFT);
threshold=0.05*max(sumFFT);
[aa,var]=sort(abs(sumFFT-threshold));
var=var(1:40); %40 big enough not to cut. ##to improve
thX1=fx(min(var));
thX2=fx(max(var));
end