Szereg Fouriera - drobna pomoc?

0

Mam prosty program do napisania, kalkuluje szereg fouriera, pakuje współczynniki w array i drukuje dane do stdout, które potem mają być wyświetlone w xmgrace. Współczynniki obliczane są metodą simpsona.

to co na razie napisałem działa dobrze, poza tym, że funkcja obliczona z fouriera ma jakby mniejszą amplitudę.
Nie mogę znaleźć błędu, byłbym bardzo wdzięczny za wszelką pomoc, wskazówki.

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
 
#define PI    3.14159265358979323846
#define PI2   6.28318530717958647692
#define PI_2  1.57079632679489661923
#define INVPI 0.31830988618379067153
 
double f1 (double x, int k)
{
        return x*sin(x*x);
}
 
double f2 (double x, int k)
{
        return f1(x, k)*cos(k*x);
}
 
double f3 (double x, int k)
{
        return f1(x, k)*sin(k*x);
}
 
/* Metoda Simpsona*/
double integrate( double (*func)(double, int), int k, int n)
        {
 
        int i;
        static double itg;
        double h;
 
        h=(PI2/n);
 
        for( i=1, itg=0 ; i<=n ; i++)
                {
                itg+=( (h/6)*(func(-PI+i*h,k) + func(-PI+(i-1)*h,k) + 4*func(0.5 * (-PI + i*h + -PI + (i-1)*h),k)) );
                }
 
        return itg;
        }
 
/* Obliczanie wspolczynnikow*/
double fourier ( double (*func1)(double, int), double (*func2)(double, int), double (*func3)(double, int), int n, int nmax, double *a, double *b)
{
        int k;
 
        b[0]=0;
        a[0]=INVPI * integrate ( func1 , 1 , n );
 
        for (k=1 ; k<=nmax ; k++)
                {
                a[k]=INVPI * integrate ( func2 , k , n );
                b[k]=INVPI * integrate ( func3 , k , n );
                }
 
        return *a, *b;
}
 
void fourier_print ( double (*func)(double, int), int nmax, double *a, double *b, int ndiv )
{                        
        int i, n;
        static double frr;
 
        for(i=0 ; i<=ndiv ; i++)
        {
                fprintf( stdout, "%2.6lf %2.6lf\n", (-PI + i*(PI2/ndiv)) , func(-PI + i*(PI2/ndiv), 0) );
        }
 
        fprintf(stdout, "&\n");
 
        for(i=0 ; i<=ndiv ; i++)
        {
                for(n=1, frr=(a[0]/2) ; n<=nmax ; n++)
                {
                        frr+=( a[n]*cos(n*(-PI + i*(PI2/ndiv))) + b[n]*sin(n*(-PI + i*(PI2/ndiv))) );
                }
                fprintf( stdout, "%2.6lf %2.6lf\n", (-PI + i*(PI2/ndiv)) , frr );
        }
}
 
 
int main ( int argc, char *argv[] )
{
 
        int nmax, ndiv, n;
        double *a;
        double *b;
 
        fprintf(stderr, "Enter n: ");            scanf("%i", &n);
        fprintf(stderr, "Enter nmax: ");       scanf("%i", &nmax);
        fprintf(stderr, "Enter ndiv: ");        scanf("%i", &ndiv);
 
        a =  malloc ( (nmax+1) * sizeof(double) );
        b =  malloc ( (nmax+1) * sizeof(double) );
 
        fourier ( f1, f2, f3, n, nmax, a, b);
 
        fourier_print ( f1 , nmax, a, b, ndiv );
 
        return EXIT_SUCCESS;
 
}
0

Mniejszą amplitudę w stosunku do czego? Być może ta metoda normalizuje FFT przez współczynnik sqrt(N), N - długość wektora. Spróbuj pomnożyć przez sqrt(N) i zobacz czy będzie to samo.

0

O ile mniejszą? Czy ta różnica zależy od n? Jeśli dla n->inf różnica dąży do zera to znaczy, że to jest po prostu błąd całkowania, którego nie da się uniknąć.
Spróbuj zobaczyć jaka jest ta różnica, gdy n jest potęgą 2.

0

Rzeczywiście, przy 'n' rzędu >10000 błąd znika, niestety docelowo muszę jeszcze wprowadzić ograniczniki dla wartości n, ndiv i nmax. Dziękuję za pomoc i fatygę.

1 użytkowników online, w tym zalogowanych: 0, gości: 1