9
|
1 |
%======================================================================
|
|
2 |
% V E C _ C O H E R E N C E . M
|
|
3 |
% doc: Wed Aug 8 09:01:59 2012
|
|
4 |
% dlm: Wed Aug 8 10:19:53 2012
|
|
5 |
% (c) 2012 A.M. Thurnherr
|
|
6 |
% uE-Info: 58 27 NIL 0 0 72 10 2 4 NIL ofnI
|
|
7 |
%======================================================================
|
|
8 |
|
|
9 |
% calculate coherence between two 2-D vector time-series or profiles
|
|
10 |
% - time or space is coded in 1st vector dimension
|
|
11 |
% - code based on [ladcp_cohere.m] by X. Liang
|
|
12 |
% - XL had hard-coded nsegs=4
|
|
13 |
|
|
14 |
function coh = vec_coherence(v1,v2,nsegs,conf)
|
|
15 |
|
|
16 |
valid = find(isfinite(v1(:,2)));
|
|
17 |
v1 = v1(find(isfinite(v1(:,2))),:);
|
|
18 |
v2 = v2(find(isfinite(v2(:,2))),:);
|
|
19 |
|
|
20 |
rez1 = median(diff(v1(:,1)));
|
|
21 |
rez2 = median(diff(v2(:,1)));
|
|
22 |
|
|
23 |
if (rez1 ~= rez2)
|
|
24 |
error('incompatible vector resolutions');
|
|
25 |
end
|
|
26 |
|
|
27 |
first = max(min(v1(:,1)),min(v2(:,1)));
|
|
28 |
last = min(max(v1(:,1)),max(v2(:,1)));
|
|
29 |
|
|
30 |
cv1 = v1(find(v1(:,1)>=first&v1(:,1)<=last),2) + sqrt(-1)*v1(find(v1(:,1)>=first&v1(:,1)<=last),3);
|
|
31 |
cv2 = v2(find(v2(:,1)>=first&v2(:,1)<=last),2) + sqrt(-1)*v2(find(v2(:,1)>=first&v2(:,1)<=last),3);
|
|
32 |
|
|
33 |
if (length(cv1) ~= length(cv2))
|
|
34 |
error('incompatible vectors');
|
|
35 |
end
|
|
36 |
|
|
37 |
NFFT = floor(length(cv1)/nsegs) ; % segment data
|
|
38 |
freq = -1/(2*rez1):1/((NFFT-1)*rez1):1/(2*rez1); % spectral frequencies
|
|
39 |
coh.conf = conf;
|
|
40 |
coh.clim = sqrt(1-(1-conf)^(1/(nsegs-1))); % Confidence limit (Data analysis methods in physical oceanography, p488)
|
|
41 |
|
|
42 |
s1 = fftshift(fft(cv1,NFFT)); % power spectra
|
|
43 |
s2 = fftshift(fft(cv2,NFFT));
|
|
44 |
|
|
45 |
s12 = conj(s1).*s2; %Inner-cross spectrum % cross spectra
|
|
46 |
s21 = fliplr(s1).*s2; %outer-cross spectrum
|
|
47 |
s11 = abs(s1).^2;
|
|
48 |
s22 = abs(s2).^2;
|
|
49 |
|
|
50 |
for i = 1:NFFT/nsegs % frequency-average spectra
|
|
51 |
as12(i) = mean(s12((i-1)*nsegs+1:i*nsegs));
|
|
52 |
as21(i) = mean(s21((i-1)*nsegs+1:i*nsegs));
|
|
53 |
as11(i) = mean(s11((i-1)*nsegs+1:i*nsegs));
|
|
54 |
as22(i) = mean(s22((i-1)*nsegs+1:i*nsegs));
|
|
55 |
coh.freq(i) = mean(freq((i-1)*nsegs+1:i*nsegs));
|
|
56 |
end
|
|
57 |
|
|
58 |
coh.period = 1 ./ coh.freq;
|
|
59 |
coh.mag = sqrt(abs(as12).^2./(as11.*as22));
|
|
60 |
coh.phase = atan2(-imag(as12),real(as12))*(180/pi);
|
|
61 |
|
|
62 |
figure
|
|
63 |
subplot(211)
|
|
64 |
plot(coh.freq,coh.mag)
|
|
65 |
hold on
|
|
66 |
line([coh.freq(1) coh.freq(end)],[coh.clim coh.clim],'color','r')
|
|
67 |
xlim([coh.freq(1) coh.freq(end)])
|
|
68 |
xlabel('Frequency/Wavenumber [cycles/unit]')
|
|
69 |
ylabel('Inner Rotary Coherence')
|
|
70 |
grid on
|
|
71 |
text(coh.freq(fix(0.8*length(coh.freq))), coh.clim+0.05, sprintf('%d%%',100*coh.conf));
|
|
72 |
|
|
73 |
subplot(212)
|
|
74 |
plot(coh.freq,coh.phase)
|
|
75 |
xlim([coh.freq(1) coh.freq(end)])
|
|
76 |
ylim([-180 180])
|
|
77 |
xlabel('Frequency/Wavenumber [cycles/unit]')
|
|
78 |
ylabel('Phase angle [deg]')
|
|
79 |
grid on
|
|
80 |
|