谱分析中的Ogive表达及其matlab实现

来源：互联网发布：网络串口服务器编辑：程序博客网时间：2024/06/05 05:23

Ogive，即谱/谐谱的累积积分的表达，如果对谱/谐谱不理解的同学，请参考傅里叶变换等相关的内容

Ogive图常用于分析特定频率对某变量方差或协方差的贡献，早期用于气象方法的湍流谱分析中，用于观测不同频率的湍流谱对湍流通量的贡献及影响，从而选择合适的频率区间

具体公式如下：

上式是分析两个变量的协谱，其中w是变量1，x是变量2，Cowx是这两个变量的傅里叶变换后的协谱，具体实现方式也给出了，Ogive的典型形状如下图：

公式下第一张图是协谱图，第二张就是Ogive curve，大致是S型，横坐标从右到左为频率重高频到低频的积分，低频部分收敛，表示对这两个变量的协方差没有多少贡献了

这里需要注意的是，Ogvie分析之前，需要对数量进行去趋势处理，detrend()，这个函数即可

function [Ogive_norm,freq_domain,Ogive,Cosp] = OgiveFunV2( X,Y,Hz,taper )%   Ogive analysis%   input:%   X-变量1%   Y - 变量2%   Hz - 频率%   taper - 窗函数%   edges at the boundaries removing the corresponding discontinuity in the%   infinitely periodic extension of the segment (Kaimal and Kristensen, 1991)%   taper - 'none';'hamming';'Bell'%   output:%   Ogive_norm: 标准化的Ogive，为啥要标准呢，当然是为了方便对比%   freq_domain: frequency domain correspond to Ogive value%   Ogive:no-normalized Ogive%   Cosp:cospectral W'X'

%%variable should detrend delete meso-flow narginchk(3,4);nargoutchk(1,5);if nargin == 3    taper = 'none';endif size(X,2) ~= 1    X = X';    Y = Y';endL = length(X); % Evaluate number of data pointsswitch taper    case 'none'        % do nothing    case 'hamming'

        %对原始数据加窗，是为了降低谱泄漏现象，不理解谱泄漏的，请参考傅里叶变换相关的内容吧或问问百度，这里给了两个选择海明窗和贝尔窗        % ref Kristoffer Aslstad's master's thesis <Applying the Eddy Covariance Method Under Difficult Conditions>        w_hamming = hamming(L,'periodic');        Xk_var = X.*w_hamming;        %the result of traper will both damp the variance and introduce a        %residual mean in the LDT time series        %first remove any residual mean introduced        Xk_var = Xk_var - mean(Xk_var);        %second calculate compensation factor recover the variance of LDT        %data        Cf = sqrt(var(X)/var(Xk_var));        Xk_var = Xk_var.*Cf;                Yk_var = Y.*w_hamming;        Yk_var = Yk_var - mean(Yk_var);        Cf = sqrt(var(Y)/var(Yk_var));        Yk_var = Yk_var.*Cf;                X = Xk_var;        Y = Yk_var;    case 'bell'        w_bell = tukeywin(L,0.25);        Xb_var = X.*w_bell;        Xb_var = Xb_var - mean(Xb_var);        %second calculate compensation factor recover the variance of LDT        %data        Cf = sqrt(var(X)/var(Xb_var));        Xb_var = Xb_var.*Cf;                Yb_var = Y.*w_bell;        Yb_var = Yb_var - mean(Yb_var);        Cf = sqrt(var(Y)/var(Yb_var));        Yb_var = Yb_var.*Cf;                X = Xb_var;        Y = Yb_var;end%plot for debugif mod(L(1),2)==0    f = Hz*linspace(0,1,L+1)';    f = f(1:end-1);else    f = Hz*linspace(0,1,L)';enddf   = f(2); % Frequency intervalfreq_domain = f(2:floor(L/2)+1); % Frequency up to Nyquist frequency

%一下两行是快速傅里叶变换y1    = fft(X)./L; % Discrete Fourier transformy2    = fft(Y)./L;hl = ones(L,1);%Cospectrum计算两个变量的协谱%sum(Co) is equal to cov(W,X,1)Co    = real(conj(hl).*conj(y1).*y2);if mod(L,2)==1    Cosp = 2.*Co(2:floor(L/2)+1,:);else    Cosp=nan(L/2,1);    Cosp(1:L/2-1,:) = 2.*Co(2:L/2,:);    Cosp(L/2,:)     = Co(L/2+1,:);endfreq_domain=real(freq_domain);               %raw natural frequencyCosp=real(Cosp);         %raw natural spectraOgive = cumsum(Cosp);%normalized Ogive%use SAM这里的SAM等同与两个变量的协方差SAM = sign(max(Ogive) + min(Ogive))*max(abs(Ogive));Ogive_norm = Ogive ./ SAM;end

0 0