dspl/group__delay_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 <string.h>

#include <math.h>

#include "dspl.h"


#ifdef DOXYGEN_ENGLISH


#endif

#ifdef DOXYGEN_RUSSIAN


#endif

int DSPL_API group_delay(double* pb, double* pa, int ord, int flag,

                         double* w, int n, double* tau)

{

    double a, b, c, d, da, db, dc, dd, f, e;

    int t, m;


    double *qa = NULL;


    if(!pb || !w || !tau || (!pa && (flag & DSPL_FLAG_ANALOG)))

        return ERROR_PTR;

    if(ord < 1)

        return ERROR_FILTER_ORD;

    if(n < 1)

        return ERROR_SIZE;


    if(pa)

      qa = pa;

    else

    {

        qa = (double*)malloc((ord+1) * sizeof(double));

        memset(qa, 0, (ord+1) * sizeof(double));

        qa[0] = 1.0;

    }


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

    {

        a = b = c = d = da = db = dc = dd = 0.0;

        if(flag & DSPL_FLAG_ANALOG)

        {

            for(m = 0; m < ord+1; m+=4)

            {

                a  += pb[m] * pow(w[t], (double)m);

                c  += qa[m] * pow(w[t], (double)m);

                da += pb[m] * (double) m * pow(w[t], (double)(m-1));

                dc += qa[m] * (double) m * pow(w[t], (double)(m-1));

            }

            for(m = 2; m < ord+1; m+=4)

            {

                a  -= pb[m] * pow(w[t], (double)m);

                c  -= qa[m] * pow(w[t], (double)m);

                da -= pb[m] * (double) m * pow(w[t], (double)(m-1));

                dc -= qa[m] * (double) m * pow(w[t], (double)(m-1));

            }


            for(m = 1; m < ord+1; m+=4)

            {

                b  += pb[m] * pow(w[t], (double)m) ;

                d  += qa[m] * pow(w[t], (double)m) ;

                db += pb[m] * (double) m * pow(w[t], (double)(m-1)) ;

                dd += qa[m] * (double) m * pow(w[t], (double)(m-1)) ;

            }


            for(m = 3; m < ord+1; m+=4)

            {

                b  -= pb[m] * pow(w[t], (double)m) ;

                d  -= qa[m] * pow(w[t], (double)m) ;

                db -= pb[m] * (double) m * pow(w[t], (double)(m-1)) ;

                dd -= qa[m] * (double) m * pow(w[t], (double)(m-1)) ;

            }


        }

        else

        {

            for(m = 0; m < ord+1; m++)

            {

                a += pb[m] * cos(w[t]*(double)m);

                b -= pb[m] * sin(w[t]*(double)m);

                c += qa[m] * cos(w[t]*(double)m);

                d -= qa[m] * sin(w[t]*(double)m);


                da -= pb[m] *(double)m * sin(w[t]*(double)m);

                db -= pb[m] *(double)m * cos(w[t]*(double)m);

                dc -= qa[m] *(double)m * sin(w[t]*(double)m);

                dd -= qa[m] *(double)m * cos(w[t]*(double)m);

            }

        }


        f = da * c + a * dc + db * d + b * dd;

        e = db * c + b * dc - da * d - a * dd;

        tau[t] = (f * (b * c - a * d) - e * (a * c + b * d)) /

                 ((a * a + b * b) * (c * c + d * d));

    }


    if(qa != pa)

      free(qa);


    return RES_OK;

}