dspl/psd__welch_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_welch(double* x, int n,

                      int win_type, double win_param,

                      int nfft, int noverlap, fft_t* pfft, double fs,

                      int flag, double* ppsd, double* pfrq)

{

    int err, pos, cnt, k;

    double *pdgr = NULL;

    double *tmp = NULL;

    double *w = NULL;

    fft_t *ptr_fft = NULL;

    double wn;


    pos = cnt = 0;


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

    if(!pdgr)

        return ERROR_MALLOC;


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

    if(!tmp)

        return ERROR_MALLOC;


    /* window malloc */

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

    if(!w)

    {

        err = ERROR_MALLOC;

        goto exit_label;

    }


    /* create window */

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

    if(err != RES_OK)

        goto exit_label;


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

    wn = 0.0;

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

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


    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)

    {

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

            tmp[k] = x[pos+k] * w[k];

        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];

        pos += noverlap;

        cnt++;

    }


    if(pos < n)

    {


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

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

            tmp[k] = x[pos+k] * w[k];


        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];


        cnt++;

    }


    /* 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)

    {

        printf("wn = %f\n", wn);

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

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

    }

    else

    {

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

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

    }


exit_label:

    if(pdgr)

        free(pdgr);

    if(tmp)

        free(tmp);

    if(w)

        free(w);

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

    {

        fft_free(ptr_fft);

        free(ptr_fft);

    }

    return err;

}