#include <stdio.h>

#define N 8 /* Nは2のべき乗である事 */

void rdft(int, int, double *, int *, double *);

void main(void)
{
    int i,j,rgsi,cy=0,radix=10000,ip[N/2+2],ans[N+2];
    double scale,a[N+2],b[N+2],c[N+2],w[N/2+2],rgs,gr;
    
    for(i=0;i<N;i++) a[i]=0.0;
    a[N-3]=9856.0;
    a[N-2]=1767.0;
    a[N-1]=3431.0;
    
    for(i=0;i<N;i++) b[i]=0.0;
    b[N-3]=7123.0;
    b[N-2]=7547.0;
    b[N-1]=3872.0;
    
    ip[0]=0;
    rdft(N, 1, a, ip, w);
    rdft(N, 1, b, ip, w);
    
    c[0]=a[0]*b[0];
    c[1]=a[1]*b[1];
    for(i=2;i<N;i+=2){
        c[i  ]=a[i]*b[i]-a[i+1]*b[i+1];
        c[i+1]=a[i+1]*b[i]+a[i]*b[i+1];
    }
    
    rdft(N,-1, c, ip, w);
    
    scale=2.0/N;
    gr=1.0/radix;
    ans[0]=0;
    for(i=N-1;i>0;i--){
        rgs=c[i-1]*scale+0.5+cy;
        cy=rgs*gr;
        rgsi=rgs;
        ans[i]=rgsi-cy*radix;
    }
    
    for(i=0;i<N;i++) printf("%04d",ans[i]);
    printf("\n");
}

戻る