dspl/psd__bartlett_8c_source.html

/*

* Copyright (c) 2015-2024 Sergey Bakhurin

* Digital Signal Processing Library [http://dsplib.org]

*

* This file is part of libdspl-2.0.

*

* is free software: you can redistribute it and/or modify

* it under the terms of the GNU Lesser    General Public License as published by

* the Free Software Foundation, either version 3 of the License, or

* (at your option) any later version.

*

* DSPL is distributed in the hope that it will be useful,

* but WITHOUT ANY WARRANTY; without even the implied warranty of

* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.    See the

* GNU General Public License for more details.

*

* You should have received a copy of the GNU Lesser General Public License

* along with Foobar.    If not, see <http://www.gnu.org/licenses/>.

*/


#include <stdlib.h>

#include <stdio.h>

#include <string.h>

#include <float.h>

#include "dspl.h"


#ifdef DOXYGEN_ENGLISH


#endif

#ifdef DOXYGEN_RUSSIAN


#endif

int DSPL_API psd_bartlett(double* x, int n, int nfft,

                          fft_t* pfft, double fs,

                          int flag, double* ppsd, double* pfrq)

{

    int err, pos, k;

    double *pdgr = NULL;

    double *tmp = NULL;

    fft_t *ptr_fft = NULL;


    pos = 0;


    pdgr = (double*)malloc(nfft * sizeof(double));

    if(!pdgr)

        return ERROR_MALLOC;


    if(!pfft)

    {

        ptr_fft = (fft_t*)malloc(sizeof(fft_t));

        memset(ptr_fft, 0, sizeof(fft_t));

    }

    else

        ptr_fft = pfft;


    memset(ppsd, 0, nfft * sizeof(double));

    while(pos + nfft <= n)

    {

        err = fft_mag(x + pos, nfft, ptr_fft, fs,

                      flag & DSPL_FLAG_FFT_SHIFT, pdgr, NULL);

        if(err != RES_OK)

            goto exit_label;

        for(k = 0; k < nfft; k++)

            ppsd[k] += pdgr[k];

        pos += nfft;

    }


    if(pos < n)

    {

        tmp = (double*)malloc(nfft * sizeof(double));

        if(!tmp)

        {

            err = ERROR_MALLOC;

            goto exit_label;

        }

        memset(tmp ,0, nfft * sizeof(double));

        memcpy(tmp, x + pos, (n - pos)*sizeof(double));


        err = fft_mag(tmp, nfft, ptr_fft, fs,

                      flag & DSPL_FLAG_FFT_SHIFT, pdgr, NULL);

        if(err != RES_OK)

            goto exit_label;


        for(k = 0; k < nfft; k++)

            ppsd[k] += pdgr[k];

    }


    /* fill frequency */

    if(pfrq)

    {

        if(flag & DSPL_FLAG_FFT_SHIFT)

            if(n%2)

                err = linspace(-fs*0.5 + fs*0.5/(double)nfft,

                                fs*0.5 - fs*0.5/(double)nfft,

                                n, DSPL_SYMMETRIC, pfrq);

            else

                err = linspace(-fs*0.5, fs*0.5, nfft, DSPL_PERIODIC, pfrq);

        else

            err = linspace(0, fs, nfft, DSPL_PERIODIC, pfrq);

    }


    /* scale magnitude */

    if(flag & DSPL_FLAG_LOGMAG)

    {

        for(k = 0; k < nfft; k++)

            ppsd[k] = 10.0 * log10(ppsd[k] / (double)n / fs);

    }

    else

    {

        for(k = 0; k < nfft; k++)

            ppsd[k] /= (double)n * fs;

    }


exit_label:

    if(pdgr)

        free(pdgr);

    if(tmp)

        free(tmp);

    if(ptr_fft && (ptr_fft != pfft))

    {

        fft_free(ptr_fft);

        free(ptr_fft);

    }

    return err;

}