L16084096 RJJ

.vscode/c_cpp_properties.json

@@ -0,0 +1,21 @@
+{
+    "configurations": [
+        {
+            "name": "Win32",
+            "includePath": [
+                "${workspaceFolder}/**"
+            ],
+            "defines": [
+                "_DEBUG",
+                "UNICODE",
+                "_UNICODE"
+            ],
+            "windowsSdkVersion": "10.0.19041.0",
+            "compilerPath": "C:/Program Files (x86)/Microsoft Visual Studio/2019/BuildTools/VC/Tools/MSVC/14.29.30133/bin/Hostx64/x64/cl.exe",
+            "cStandard": "c17",
+            "cppStandard": "c++17",
+            "intelliSenseMode": "windows-msvc-x64"
+        }
+    ],
+    "version": 4
+}

5_practice.c

@@ -9,32 +9,43 @@ int main() {
 	time_t t;
 
 	// Input the number N
-
+	printf("Enter N = ");
+	scanf("%d",&N);
 
 	// Locate the memory for x, yr, yi;
-
-	// Initial setting for x, for example, x[k] = k
-
-
+	x = (double *) malloc(N*sizeof(double));
+	yr = (double *) malloc(N*sizeof(double));
+	yi = (double *) malloc(N*sizeof(double));
 
+	// Initial setting for x, for example, x[k] = k
+	for(int i = 0; i<N; i++){
+		x[i] = i;
+		yr[i] = 0.0;
+		yi[i] = 0.0;
+	}
 
 	t = clock();	
 	// yr[n]+i*yi[n] = sum(exp(-i 2 Pi k n / N)*x[k], k=0..N-1), n=0..N-1 
-
-
-
-
-
-
+	//					cos(2 PI k n / N)		i sin(2 PI k n / N)
+	double tmp = 2*M_PI/N;
+	double theta;
+	for (int i = 1; i<N; i++){
+ 		for(int j = 1; j<N; j++){
+			theta = tmp*i*j;
+			yr[i] += cos(theta)*x[i];
+			yi[i] -= sin(theta)*x[i];
+		}
+	}
 
 	// output the results
 	t = clock() - t;
 	printf("%d ms for discrete Fourier Transform of %d elements\n", t, N);
 
 	// free the memory located by x, yr, yi
-
-
-
-	return 100;
+	free(yi);
+	free(yr);
+	free(x);
+
+	return 0;
 }
 

FFTP.c

@@ -0,0 +1,130 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include "complex.h"
+
+int roundup(int a, int b){
+    int c;
+    if(a%b){
+        a += (b-a%b);
+    }
+    c = a/b;
+    return c;
+}
+
+complex FFTP(complex *y, complex *x, int N, int p){
+
+    //declare some variables 
+    int i, j, k, Np = N/p;     //i,j,k for index using
+
+    //prime factorization
+    //delcare some variables
+    int nextp, nextN, psize;
+    int *prime, *ptime;     //record prime and square number
+    int Ntmp = N, index = 0, t, ctrl, final;
+
+    //memory allocation for prime and its square number
+    psize = roundup(N,2);
+    prime = (int *) malloc( psize * sizeof(int));
+    ptime = (int *) malloc( psize * sizeof(int));
+
+    for(i = 2; i < psize; i++){
+        t = 0;
+        while(Ntmp%i == 0){
+            Ntmp /= i;
+            t++;                //sum up square number
+            ctrl = 1;           //switch the position of the array
+        }
+
+        if(ctrl == 1){
+            prime[index] = i;   //prime        
+            ptime[index] = t;   //square number of the prime resp.
+            index++;
+        }
+        ctrl = 0;
+    }
+    if(!index){
+        nextp = N;
+        final = 1;              //as N = p, do the final FFT
+        prime[index] = N;
+        ptime[index] = 1;
+    }else if(ptime[index-1] > 1){
+        nextp = prime[index-1];
+    }else{
+        nextp = prime[index-2];
+    }
+
+    //FFTP
+    if(Np == 1 || final){
+        p = N;
+        complex *w;
+        w = (complex *) malloc( p * sizeof(complex) );
+        for(i = 0; i < p; i++){
+            w[i] = wtrans(N,i);
+        }
+        for(i = 0; i < p; i++){
+            y[i].Re = 0;
+            y[i].Im = 0;
+            for(j = 0; j < p; j++){
+                y[i] = Cadd(y[i], Cmulti( x[j] ,w[(i*j)%N] ) );
+            } 
+        }
+        free(w);
+    }else{
+
+        //memory allocation for subseq
+        nextN = N / nextp;
+        complex **xsub, **ysub;
+        xsub = (complex **) malloc( nextp * sizeof(complex) );
+        ysub = (complex **) malloc( nextp * sizeof(complex) );
+        for(i = 0; i < nextp; i++){
+            xsub[i] = (complex *) malloc( nextN * sizeof(complex) );
+            ysub[i] = (complex *) malloc( nextN * sizeof(complex) );
+        }
+
+        //assign x subseq and y
+        for(i = 0; i < nextN; i++){
+            for(j = 0; j < nextp; j++){
+                xsub[i][j] = x[nextp*j + i];
+                y[i*nextN + j].Re = 0;
+                y[i*nextN + j].Im = 0;
+            }
+        }
+
+        //do the subseq FFT
+        for(i = 0;i < nextp; i++){
+            FFTP( *(ysub+i), *(xsub+i), nextN, nextp);
+        }
+
+        //combine the subseq 
+        complex *w;
+        w = (complex *) malloc( nextp * sizeof(complex));
+        for( j = 0; j < nextp; j++){
+            w[j] = wtrans(N,i);
+        }
+        for(i = 0; i < nextp; i++){
+            for( j = 0; j < nextN; j++){
+                for(k = 0; k < nextp; k++){
+                    y[ nextp * i + j ] = Cadd( y[ nextp * i + j ] , \
+					Cmulti(ysub[k][j], w[ ( k* (i*N/p+j) ) %N] ) );
+                }
+            }
+        }
+        //free every memory
+        free(w);
+        for(i = 0; i < nextp; i++){
+            free(xsub[i]);
+            free(ysub[i]);
+        }
+        free(xsub);free(ysub);
+
+    }
+
+    free(prime);free(ptime);
+
+    //end
+    complex c;
+    c.Re = 0;
+    c.Im = 0;
+    return c;
+}

FFT_235.c

@@ -0,0 +1,103 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <time.h>
+#include "fft.h"
+
+//declare the function FFT
+int FFT(double *yr, double *yi, double *xr, double *xi, int N);
+
+int main(){
+
+    //declare variables
+    //Suppose N is the total number, pqr are square of 235 respec.
+    //Suppose 
+    int k, n, N, p, q, r, ctrl;
+    double *xr, *xi, *yr, *yi, t, a, b, an, bn, c, s;
+    time_t T;
+
+    xr = (double *) malloc(N*sizeof(double));
+    xi = (double *) malloc(N*sizeof(double));
+    yr = (double *) malloc(N*sizeof(double));
+    yi = (double *) malloc(N*sizeof(double));
+
+    printf("Memory Done! \n");
+
+    for(k = 0;k<5;++k){
+        printf("%d: %f \n",k,xr[k]);
+    }
+
+    T = clock();
+
+    a = cos(2 * M_PI / N); b = sin(2 * M_PI / N);
+    an = 1; bn = 0;
+    for(n = 0; n < N; n++){
+        yr[n] = 0.0;
+        yi[n] = 0.0;
+        c = 1; s = 0;
+        for(k = 0; k<N; ++k){
+            yr[n] += c*xr[k];
+            yi[n] -= s*xr[k];
+            t = c * an - bn * b;
+            s = c * bn + s * an;
+            c = t;
+        }
+        t = an * a - bn * b;
+        bn = an * b + bn * a;
+        an = t;
+    }
+
+    T = clock() - T;
+    printf("%d ms for DFT of %d elements\n",T,N);
+
+    free(xr);free(xi);free(yr);free(yi);
+
+    return 0;
+}
+
+int FFT(double *yr, double *yi, double *xr, double *xi, int N){
+    if(N == 2){
+        yr[0] = xr[0]+xr[1];
+        yi[0] = xi[0]+xi[1];
+        yr[1] = xr[0]-xr[1];
+        yr[1] = xi[0]-xi[1];
+    } else {
+        int k;
+        double *yEr, *yEi, *yOr, *yOi;
+        double *xEr, *xEi, *xOr, *xOi;
+        double c, s;
+
+        xEr = (double *) malloc((N/2)*sizeof(double));
+        xEi = (double *) malloc((N/2)*sizeof(double));
+        yEr = (double *) malloc((N/2)*sizeof(double));
+        yEi = (double *) malloc((N/2)*sizeof(double));
+        xOr = (double *) malloc((N/2)*sizeof(double));
+        xOi = (double *) malloc((N/2)*sizeof(double));
+        yEr = (double *) malloc((N/2)*sizeof(double));
+        yOi = (double *) malloc((N/2)*sizeof(double));
+
+        for(k = 0; k<N/2; ++k){
+            xEr[k] = xr[2*k];
+            xEi[k] = xi[2*k];
+            xOr[k] = xr[2*k+1];
+            xOi[k] = xi[2*k+1];
+
+        }
+        FFT(yEr, yEi, xEr, xEi, N/2);
+        FFT(yOr, yOi, xOr, xOi, N/2);
+        for(k = 0; k < N/2; ++k){
+            c = cos(2*M_PI*k/N);
+            s = -sin(2*M_PI*k/N);
+            yr[k] = yEr[k] + (yOr[k]*c - yOi[k]*s);
+            yi[k] = yEi[k] + (yOr[k]*s + yOi[k]*c);
+            yr[N/2+k] = yEr[k] - (yOr[k]*c - yOi[k]*s);
+            yi[N/2+k] = yEi[k] - (yOr[k]*s + yOi[k]*c);
+        }
+
+        free(xEr);free(xEi);
+        free(xOr);free(xOi);
+        free(yEr);free(yEi);
+        free(yOr);free(yOi);
+    }
+    return 0;
+}

Readme_FFT.txt

@@ -0,0 +1,3 @@
+�ڪ�FFT�]�p���Q�k�O�����N��N�A������]�Ƥ��ѡA�̾ڳ̤j���p_max���l�ǦC���ΡA�ҥH�N�|��p_max�Ӫ��׬�N/p_max���l�ǦC�A���o��o�Ǥl�ǦC�A���@�˪��B�J�A����C�Ӥl�ǦC�������ƶq����ƫK�}�l���̧C�h���ť߸��ۥ[�A���۩��W���|�^�h�o��̫᪺���סC
+
+�{�b�ɶ��O4/29 23:55�A�ثe���i�׬O���]�A���O�]�X�Ӽƾڤ��ӹ�A�ثe�����O�bFFT�|�N���ɭԥX�F���D�C