LPSD 功率谱估计

% Use LPSD mothod to calculate power spectral density
% Usage:
%    h = iLPSD(data,fs,_)
%    [pxx,f] = iLPSD(data,fs,_)
% Required paramters:
%    data --- input data, processed by column
%    fs   --- sample frequency, unit: Hz
% Optional paramters:
%    Jdes --- desired frequency points, default: 1000
%    Kdes --- desired number of segments, default: 100
%    xi   --- overlap ratio, default: 0.5
%    win  --- window function handle, default: @hann
%    par  --- enable parallel compution, default: false
%    type --- type of the sepctrum: PSD|RMS|Amp, default: PSD
%    Extra paramters to control plotting are allowed (e.g. 'LineWidth')
% Output paramters
%    h    --- line handle of the plotted PSD
%    pxx  --- [PSD] one-sided PSD, unit: *^2/Hz
%             [RMS| Amp] RMS or amplitude, unit: *
%    f    --- frequency points related to PSD points, unit: Hz
%
% Ref: Improved spectrum estimation from digitized time series on a logarithmic frequency axis
%      Article DOI: 10.1016/j.measurement.2005.10.010
%
% For more information, see <a href="https://ichunyu.github.io/helps/functions/ilpsd"
% >online documentation</a>.

% XiaoCY 2019-02-23, first version (date from PPT report)
% XiaoCY 2020-04-21, support matrix data: treated as columns
% XiaoCY 2021-11-06, add inputParser to parse parameters
% XiaoCY 2022-03-28, add parameters: 'parallel', 'type'
% XiaoCY 2022-05-08, bug fix: window function in parallel mode

%% Main
function varargout = iLPSD(varargin)

    % parse input paramters
    p = inputParser;
    p.FunctionName = 'iLPSD';
    p.KeepUnmatched = true;

    p.addRequired('data');
    p.addRequired('fs');
    p.addParameter('Jdes',1000);
    p.addParameter('Kdes',100);
    p.addParameter('xi',0.5);
    p.addParameter('window',@hann);
    p.addParameter('parallel',false);
    p.addParameter('type','PSD');
    p.parse(varargin{:});

    % unpack necessary parameters
    data = p.Results.data;
    fs = p.Results.fs;
    Jdes = p.Results.Jdes;
    Kdes = p.Results.Kdes;
    xi = p.Results.xi;
    win = p.Results.window;
    type = validatestring(p.Results.type,{'PSD','RMS','Amp'});
    opts = p.Unmatched;
    
    % ensure input data as column vector(s)
    [N,nCol] = size(data);
    if N==1 && nCol~=1
        data = data';
        N = nCol;
        nCol = 1;
    end

    [f,L,m] = getFreqs(N,fs,Jdes,Kdes,xi);

    J = length(f);
    P = zeros(J,nCol);
    if p.Results.parallel
        parfor j = 1:J
            Dj = floor((1-xi)*L(j));
            Kj = floor((N-L(j))/Dj+1);
            w = win(L(j));
            l = (0:L(j)-1)';
            W1 = cos(-2*pi*m(j)/L(j).*l);
            W2 = sin(-2*pi*m(j)/L(j).*l);
            A = zeros(1,nCol);
            for k = 0:Kj-1
                G = data(k*Dj+1:k*Dj+L(j),:);
                G = G-mean(G);        % G = detrend(G);
                G = G.*w;
                A = A + sum(G.*W1).^2+sum(G.*W2).^2;
            end
            switch type
                case 'PSD'
                    P(j,:) = A/Kj*(2/fs/sum(w.^2));
                case 'RMS'
                    P(j,:) = A/Kj*(2/sum(w)^2);
                otherwise
                    P(j,:) = 2*A/Kj*(2/sum(w)^2);
            end
        end
    else
        for j = 1:J
            Dj = floor((1-xi)*L(j));
            Kj = floor((N-L(j))/Dj+1);
            w = win(L(j));
            l = (0:L(j)-1)';
            W1 = cos(-2*pi*m(j)/L(j).*l);
            W2 = sin(-2*pi*m(j)/L(j).*l);
            A = zeros(1,nCol);
            for k = 0:Kj-1
                G = data(k*Dj+1:k*Dj+L(j),:);
                G = G-mean(G);        % G = detrend(G);
                G = G.*w;
                A = A + sum(G.*W1).^2+sum(G.*W2).^2;
            end
            switch type
                case 'PSD'
                    P(j,:) = A/Kj*(2/fs/sum(w.^2));
                case 'RMS'
                    P(j,:) = A/Kj*(2/sum(w)^2);
                otherwise
                    P(j,:) = 2*A/Kj*(2/sum(w)^2);
            end
        end
    end

    switch nargout
        case 0
            loglog(f,sqrt(P),opts);
            grid on
            xlabel('Frequency [Hz]','interpreter','latex')
            switch type
                case 'PSD'
                    ylabel('PSD [$\mathrm{*/\sqrt{Hz}}$]','interpreter','latex')
                case 'RMS'
                    ylabel('RMS [$\mathrm{*}$]','interpreter','latex')
                otherwise
                    ylabel('Amp [$\mathrm{*}$]','interpreter','latex')
            end
        case 1
            varargout{1} = loglog(f,sqrt(P),opts);
            grid on
            xlabel('Frequency [Hz]','interpreter','latex')
            switch type
                case 'PSD'
                    ylabel('PSD [$\mathrm{*/\sqrt{Hz}}$]','interpreter','latex')
                case 'RMS'
                    ylabel('RMS [$\mathrm{*}$]','interpreter','latex')
                otherwise
                    ylabel('Amp [$\mathrm{*}$]','interpreter','latex')
            end
        case 2
            if strcmp(type,'PSD')
                varargout{1} = P;
            else
                varargout{1} = sqrt(P);
            end
            varargout{2} = f;
        otherwise
            % Do Nothing
    end
end

%% Subfunctions
% get logarithmic frequency points
function [f,L,m] = getFreqs(N,fs,Jdes,Kdes,xi)
    fmin = fs/N;
    fmax = fs/2;
    r_avg = fs/N*(1+(1-xi)*(Kdes-1));

    g = (N/2)^(1/(Jdes-1))-1;

    f = zeros(Jdes,1)-1;
    L = f;
    m = f;
    j = 1;
    fj = fmin;
    while fj < fmax
        rj = fj*g;
        if rj < r_avg
            rj = sqrt(rj*r_avg);
        end
        if rj < fmin
            rj = fmin;
        end

        Lj = floor(fs/rj);
        rj = fs/Lj;
        mj = fj/rj;

        f(j) = fj;
        L(j) = Lj;
        m(j) = mj;

        fj = fj+rj;
        j = j+1;
    end
    f(f<0) = [];
    L(L<0) = [];
    m(m<0) = [];
end