Added an averaging routine after Murphy et al 2014

WaveletTransform.m

@@ -14,9 +14,10 @@
         u_wind_component;
         v_wind_component;
         power_surface;
-        u_wavelet;
-        v_wavelet;
+        uWavelet;
+        vWavelet;
         waveletScales;
+        fourierPeriod;
         coi;
     end
 
@@ -32,30 +33,36 @@
             obj.s0 = s0;
             obj.dj = dj;
             obj.dt = dt;
-
-            [obj.u_wavelet, ~, obj.waveletScales, ~] = wavelet(u_wind_component, dt, pad, dj, s0); % wave, period, scale, COI
-            [obj.v_wavelet, ~, ~, obj.coi] = wavelet(v_wind_component, dt, pad, dj, s0);
-            obj.power_surface = abs(obj.u_wavelet).^2 + abs(obj.v_wavelet).^2;
+            [obj.uWavelet, ~, obj.waveletScales, ~] = wavelet(u_wind_component, dt, pad, dj, s0); % wave, period, scale, COI
+            [obj.vWavelet, ~, ~, obj.coi] = wavelet(v_wind_component, dt, pad, dj, s0);
+            obj.power_surface = abs(obj.uWavelet).^2 + abs(obj.vWavelet).^2;
+            obj.fourierPeriod = 1.03 * obj.waveletScales;
         end
 
         function [u_wind_inverted, v_wind_inverted, v_wind_hilbert_transformed] = invertWindowedTransform(obj, windowedWaveletTransform)
-            %METHOD1 Summary of this method goes here
-            %   Detailed explanation goes here
             scale_index_1 = windowedWaveletTransform.scale_index_1;
             scale_index_2 = windowedWaveletTransform.scale_index_2;
             alt_index_1 = windowedWaveletTransform.alt_index_1;
             alt_index_2 = windowedWaveletTransform.alt_index_2;
             constant_coef = obj.dj * sqrt(obj.dt) / (0.776*pi^(-1/4)); % Magic from T&C.
             windowed_scales = obj.waveletScales(scale_index_1:scale_index_2);
-            windowed_u_wavelet = obj.u_wavelet(scale_index_1:scale_index_2, alt_index_1:alt_index_2);
-            u_wind_reconstructed = constant_coef*sum(windowed_u_wavelet ./ sqrt(windowed_scales)');
-            windowed_v_wavelet = obj.v_wavelet(scale_index_1:scale_index_2, alt_index_1:alt_index_2);
-            v_wind_reconstructed = constant_coef*sum(windowed_v_wavelet ./ sqrt(windowed_scales)');
+            windowed_u_wavelet = obj.uWavelet(scale_index_1:scale_index_2, alt_index_1:alt_index_2);
+            u_wind_reconstructed = constant_coef*sum(windowed_u_wavelet ./ sqrt(windowed_scales)', 1);
+            windowed_v_wavelet = obj.vWavelet(scale_index_1:scale_index_2, alt_index_1:alt_index_2);
+            v_wind_reconstructed = constant_coef*sum(windowed_v_wavelet ./ sqrt(windowed_scales)', 1);
             v_wind_inverted = real(v_wind_reconstructed);
             u_wind_inverted = real(u_wind_reconstructed);
             v_wind_hilbert_transformed = complex(v_wind_reconstructed);            
         end
 
+        function [uWindInverted, vWindInverted] = invertWaveletTransform(obj)
+           constantCoef = obj.dj * sqrt(obj.dt) / (0.776*pi^(1/4)); % Magic from T&C.
+           uWindInverted = constantCoef*sum(real(obj.uWavelet) ./ sqrt(obj.waveletScales)', 1); % sum over scales
+           vWindInverted = constantCoef*sum(real(obj.vWavelet) ./ sqrt(obj.waveletScales)', 1); 
+
+        end
+
+
         function [a, b, c, d] = clipWindowedTransformToValue(obj, windowedWaveletTransform, localMaxRow, localMaxCol)
             % clipWindowedTransformToValue either clips the windowed
             % transform to ``value'', or clips it to the next inflection

averageToAltitudeResolution.m

@@ -0,0 +1,32 @@
+function [averagedArray, averagedAlt] = averageToAltitudeResolution(array, altitude, resolution)
+%UNTITLED Summary of this function goes here
+%   Averages a time series of data to a given altitude resolution. Steps:
+%   walk through the altitude array, keeping a running sum of the values in
+%   "array" at the same indices. When altBegin - altEnd > resolution,
+%   divide the sum by the number of indices between altBegin and altEnd,
+%   and store it in a new array - averagedArray.
+
+averagedArray = zeros('like', array);
+averagedAlt = zeros('like', altitude);
+altBegin = altitude(1);
+altEnd = altitude(1);
+k = 0;
+s = 0;
+nPoints = 0;
+for i=1:size(altitude, 1)
+    altEnd = altitude(i);
+    s = s + array(i);
+    nPoints = nPoints + 1;
+    if (altEnd - altBegin) >= resolution
+        k = k + 1;
+        averagedArray(k) = s / nPoints;
+        averagedAlt(k) = altBegin + (altEnd - altBegin) / 2;
+        s = 0;
+        altBegin = altEnd;
+        nPoints = 0;
+    end
+end
+averagedArray = averagedArray(1:k);
+averagedAlt = averagedAlt(1:k);
+end
+

estimateParametersFromWavePacket.m

@@ -10,6 +10,7 @@
 Q = 2*mean(u.*vWavePacketHilbertTransformed);
 degreeOfPolarization = sqrt((P^2 + Q^2 + D^2)) / I;
 if Q < 0.05 || P < 0.05 || degreeOfPolarization < 0.5 || degreeOfPolarization > 1
+%if degreeOfPolarization < 0.5 || degreeOfPolarization > 1
    theta = 0;
    axialRatio = 0;
    degreeOfPolarization = 0;