1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64
|
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <fftw3.h>
#include "ran.h"
main()
{
int n;
int k1,k2;
int i,j;
double sumnoise,sumnoisek;
double *noise;
fftw_complex *noisek;
fftw_plan p1;
n=256;
noisek=fftw_malloc(n*(n/2+1)*sizeof(fftw_complex));
noise=fftw_malloc(n*n*sizeof(double));
p1=fftw_plan_dft_c2r_2d(n,n,noisek,noise,FFTW_ESTIMATE);
seed3(95539);
sumnoisek=0.0;
for(k1=0;k1<n;k1++)
{
for(k2=0;k2<n/2+1;k2++)
{
noisek[k2+k1*(n/2+1)][0]=ran3();
noisek[k2+k1*(n/2+1)][1]=ran3();
if(k2==0 || k2==n/2)
sumnoisek+=noisek[k2+k1*(n/2+1)][0]*noisek[k2+k1*(n/2+1)][0]
+noisek[k2+k1*(n/2+1)][1]*noisek[k2+k1*(n/2+1)][1];
else
sumnoisek+=2.0*(noisek[k2+k1*(n/2+1)][0]*noisek[k2+k1*(n/2+1)][0]
+noisek[k2+k1*(n/2+1)][1]*noisek[k2+k1*(n/2+1)][1]);
}
}
fftw_execute(p1);
printf("%f \n",sumnoisek);
sumnoise=0.0;
for(i=0;i<n;i++)
for(j=0;j<n;j++)
sumnoise+=noise[j+i*n]*noise[j+i*n];
printf("%f \n",sumnoise/n/n);
fftw_destroy_plan(p1);
fftw_free(noisek);
fftw_free(noise);
} |
Partager