|
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 |