# HG changeset patch # User A.M. Thurnherr # Date 1385033033 0 # Node ID 79a432124b5751891cf966be25acd77bc0056ebf # Parent dcc97ece4503a892f6015b1386bdd44f58820074 . diff --git a/LDEO_LADCP2ANTS.m b/LDEO_LADCP2ANTS.m --- a/LDEO_LADCP2ANTS.m +++ b/LDEO_LADCP2ANTS.m @@ -1,14 +1,14 @@ %====================================================================== % L D E O _ L A D C P 2 A N T S . M % doc: Sun Jan 22 15:19:00 2006 -% dlm: Mon Oct 12 22:47:23 2009 +% dlm: Mon Jun 24 10:37:03 2013 % (c) 2006 A.M. Thurnherr -% uE-Info: 155 47 NIL 0 0 72 2 2 4 NIL ofnI +% uE-Info: 30 40 NIL 0 0 72 2 2 4 NIL ofnI %====================================================================== % % export LDEO LADCP output to ANTS file % -% USAGE: LDEO_LADCP2ANTS(dr,f,p,outBaseName) +% USAGE: LDEO_LADCP2ANTS(dr,f,p,ps,outBaseName) % % HISTORY: @@ -22,8 +22,15 @@ % Jul 17, 2008: - added cruise, software, magdecl, procdir info % Apr 23, 2009: - added globarl var EXPORT_CTD_DATA % Oct 12, 2009: - adapted to new struct2ANTS +% Mar 18, 2013: - added support for global STRUCT2ANTS.verb +% Jun 24, 2013: - added blen, nbin, blnk, dist to output (DL/UL separately); V10 +% - added %depth_resolution %ADCP_superens_dz to output, requiring ps as +% additional input -function [] = LDEO_LADCP2ANTS(dr,f,p,obn) +function [] = LDEO_LADCP2ANTS(dr,f,p,ps,obn) + + global STRUCT2ANTS; % suppress diagnostic messages + STRUCT2ANTS.verb = 0; %---------------------------------------------------------------------- % INVERSE SOLUTION @@ -34,6 +41,21 @@ prof.software = p.software; prof.magdecl = p.drot; prof.procdir = pwd; + + prof.DL_bin_length = p.blen_d; + prof.DL_bins = p.nbin_d; + prof.DL_blanking = p.blnk_d; + prof.DL_bin1_dist = p.dist_d; + + if isfield(p,'nbin_u') + prof.UL_bin_length = p.blen_u; + prof.UL_bins = p.nbin_u; + prof.UL_blanking = p.blnk_u; + prof.UL_bin1_dist = p.dist_u; + end + + prof.ADCP_superens_dz = p.avdz; + prof.depth_resolution = ps.dz; prof.start_date = sprintf('%d/%02d/%02d',p.time_start(1),p.time_start(2),p.time_start(3)); prof.start_time = sprintf('%02d:%02d:%02d',p.time_start(4),p.time_start(5),p.time_start(6)); diff --git a/struct2ANTS.m b/struct2ANTS.m --- a/struct2ANTS.m +++ b/struct2ANTS.m @@ -1,9 +1,9 @@ %====================================================================== % S T R U C T 2 A N T S . M % doc: Thu Oct 20 11:48:17 2005 -% dlm: Tue Feb 21 14:30:01 2012 +% dlm: Mon Mar 18 16:17:52 2013 % (c) 2005 A.M. Thurnherr -% uE-Info: 93 0 NIL 0 0 72 2 2 4 NIL ofnI +% uE-Info: 96 20 NIL 0 0 72 2 2 4 NIL ofnI %====================================================================== % % export Matlab structure to ANTS file @@ -15,6 +15,7 @@ % - scalar strings & numbers become %PARAMs % - row and column vectors (which can be mixed) become FIELDs % - incompatible vectors, as well as other data types, are skipped +% - set global STRUCT2ANTS.verb to true to suppress diagnostic messages % HISTORY: % Oct 20, 2005: - created @@ -37,6 +38,7 @@ % Feb 21, 2012: - removed double quoting of % and $ % - manually merged two versions % - re-added diagnostic messages about skipped incompatible vectors +% Mar 18, 2013: - added global STRUCT2ANTS.verb function [] = struct2ANTS(struct,deps,ofn) @@ -64,6 +66,11 @@ function [ldef,dta] = parseStruct(struct,ldef,fpref,dta,ofn); + global STRUCT2ANTS; + if ~isfield(STRUCT2ANTS,'verb') + STRUCT2ANTS.verb = 1; + end + fname = fieldnames(struct); for i=1:length(fname) fns = char(fname(i)); @@ -86,7 +93,7 @@ elseif c==1 && length(f)==ndta ldef = sprintf('%s%s%s= ',ldef,fpref,fns); dta = [dta,f]; - else + elseif STRUCT2ANTS.verb > 0 disp(sprintf('%s: incompatible %d x %d array %s%s skipped',ofn,r,c,fpref,fns)); end end diff --git a/vec_coherence.m b/vec_coherence.m new file mode 100644 --- /dev/null +++ b/vec_coherence.m @@ -0,0 +1,80 @@ +%====================================================================== +% V E C _ C O H E R E N C E . M +% doc: Wed Aug 8 09:01:59 2012 +% dlm: Wed Aug 8 10:19:53 2012 +% (c) 2012 A.M. Thurnherr +% uE-Info: 58 27 NIL 0 0 72 10 2 4 NIL ofnI +%====================================================================== + +% calculate coherence between two 2-D vector time-series or profiles +% - time or space is coded in 1st vector dimension +% - code based on [ladcp_cohere.m] by X. Liang +% - XL had hard-coded nsegs=4 + +function coh = vec_coherence(v1,v2,nsegs,conf) + +valid = find(isfinite(v1(:,2))); +v1 = v1(find(isfinite(v1(:,2))),:); +v2 = v2(find(isfinite(v2(:,2))),:); + +rez1 = median(diff(v1(:,1))); +rez2 = median(diff(v2(:,1))); + +if (rez1 ~= rez2) + error('incompatible vector resolutions'); +end + +first = max(min(v1(:,1)),min(v2(:,1))); +last = min(max(v1(:,1)),max(v2(:,1))); + +cv1 = v1(find(v1(:,1)>=first&v1(:,1)<=last),2) + sqrt(-1)*v1(find(v1(:,1)>=first&v1(:,1)<=last),3); +cv2 = v2(find(v2(:,1)>=first&v2(:,1)<=last),2) + sqrt(-1)*v2(find(v2(:,1)>=first&v2(:,1)<=last),3); + +if (length(cv1) ~= length(cv2)) + error('incompatible vectors'); +end + +NFFT = floor(length(cv1)/nsegs) ; % segment data +freq = -1/(2*rez1):1/((NFFT-1)*rez1):1/(2*rez1); % spectral frequencies +coh.conf = conf; +coh.clim = sqrt(1-(1-conf)^(1/(nsegs-1))); % Confidence limit (Data analysis methods in physical oceanography, p488) + +s1 = fftshift(fft(cv1,NFFT)); % power spectra +s2 = fftshift(fft(cv2,NFFT)); + +s12 = conj(s1).*s2; %Inner-cross spectrum % cross spectra +s21 = fliplr(s1).*s2; %outer-cross spectrum +s11 = abs(s1).^2; +s22 = abs(s2).^2; + +for i = 1:NFFT/nsegs % frequency-average spectra + as12(i) = mean(s12((i-1)*nsegs+1:i*nsegs)); + as21(i) = mean(s21((i-1)*nsegs+1:i*nsegs)); + as11(i) = mean(s11((i-1)*nsegs+1:i*nsegs)); + as22(i) = mean(s22((i-1)*nsegs+1:i*nsegs)); + coh.freq(i) = mean(freq((i-1)*nsegs+1:i*nsegs)); +end + +coh.period = 1 ./ coh.freq; +coh.mag = sqrt(abs(as12).^2./(as11.*as22)); +coh.phase = atan2(-imag(as12),real(as12))*(180/pi); + +figure +subplot(211) +plot(coh.freq,coh.mag) +hold on +line([coh.freq(1) coh.freq(end)],[coh.clim coh.clim],'color','r') +xlim([coh.freq(1) coh.freq(end)]) +xlabel('Frequency/Wavenumber [cycles/unit]') +ylabel('Inner Rotary Coherence') +grid on +text(coh.freq(fix(0.8*length(coh.freq))), coh.clim+0.05, sprintf('%d%%',100*coh.conf)); + +subplot(212) +plot(coh.freq,coh.phase) +xlim([coh.freq(1) coh.freq(end)]) +ylim([-180 180]) +xlabel('Frequency/Wavenumber [cycles/unit]') +ylabel('Phase angle [deg]') +grid on +