dspl/psd__periodogram__cmplx_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_periodogram_cmplx(complex_t* x, int n,

                                   int win_type, double win_param,

                                   fft_t* pfft, double fs,

                                   int flag, double* ppsd, double* pfrq)

{

    double *w = NULL;

    complex_t *s = NULL;

    double u, wn;

    int err, k;

    fft_t *ptr_fft = NULL;


    if(!x || !ppsd)

        return  ERROR_PTR;


    if(n<1 )

        return ERROR_SIZE;


    if(fs < 0.0)

        return ERROR_FS;


    if(!pfft)

    {

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

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

    }

    else

        ptr_fft = pfft;


    if(win_type != DSPL_WIN_RECT)

    {

        /* Modified periodogram calculation */


        /* window malloc */

        w = (double*)malloc(n*sizeof(double));

        if(!w)

        {

            err = ERROR_MALLOC;

            goto exit_label;

        }


        /* create window */

        err = window(w, n, win_type, win_param);

        if(err != RES_OK)

            goto exit_label;


        /* window normalization wn = sum(w.^2) */

        wn = 0;

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

            wn += w[k]*w[k];


        /* signal buffer malloc */

        s = (complex_t*)malloc(n*sizeof(complex_t));

        if(!s)

        {

            err = ERROR_MALLOC;

            goto exit_label;

        }


        /* windowing */

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

        {

            RE(s[k]) = RE(x[k]) * w[k];

            IM(s[k]) = IM(x[k]) * w[k];

        }

    }

    else

    {

        /* classic periodogram without windowing */

        s = x;

        wn = (double)n;

    }


    /* calculate FFT */

    err = fft_mag_cmplx(s, n, ptr_fft, fs, flag, ppsd, pfrq);

    if(err != RES_OK)

        goto exit_label;


    if(flag & DSPL_FLAG_LOGMAG)

    {

        /* normalization in log scale */

        u = 10.0 * log10(wn * fs);

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

            ppsd[k] -= u;

    }

    else

    {

        /* normalization in linear scale */

        u = 1.0 / (wn * fs);

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

            ppsd[k] *= u;

    }


exit_label:

    if(w)

        free(w);

    if(s && s != x)

        free(s);

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

    {

        fft_free(ptr_fft);

        free(ptr_fft);

    }

    return err;

}