Finish all assignments

5_moving_average.c

@@ -2,43 +2,87 @@
 #include <stdlib.h>	// for memory allocation
 #include <time.h>	// for time calculation
 #include <math.h>	// for sine and cosine functions
-int main() {
+#include <string.h>
+
+// #define TEST
+
+#ifdef TEST
+#include <assert.h>
+#endif
+
+int main(int argc, const char *argv[]) {
 	// Declare all the variables
-	int k, n, M, N;
+	int M, N;
 	double *x, *y;
 	time_t t;
 
 	// Input the number N, M
-	printf("N M= ");	
-	scanf("%d %d", &N, &M);
-
+    printf("N = ");
+    scanf("%d", &N);	
+    printf("M = ");
+    scanf("%d", &M);
+
+    double size = 2*M + 1;
 	// Locate the memory for x, y;
-	x = (double *) malloc(N*sizeof(double));
-	y = (double *) malloc(N*sizeof(double));
+    x = (double*)malloc(sizeof(double)*N);
+    y = (double*)malloc(sizeof(double)*N);
+
+
+#ifdef TEST
+    double *y_test = (double*)malloc(sizeof(double)*N);
+    memset(y, 0, sizeof(double)*N);
+#endif
+
 	// Initial setting for x[k] = sin(2*Pi*k/N), k = 0..N-1
-	for(k=0;k<N;++k) {
-		x[k] = sin(2*M_PI*k/N);
-	}
+    for(int i = 0; i < N; ++i){
+        x[i] = sin(2*M_PI*i/(double)N);
+    }
+    memset(y, 0, sizeof(double)*N);
+
+#ifdef TEST
+    // use O(NM) moving average to perform test
+    for(int n = 0; n < N; ++n){
+        for(int k = n-M, upper_bound = n+M; k <= upper_bound; ++k){
+            if(k >= 0 && k < N){
+                y_test[n] += x[k];
+            }
+        }
+        y_test[n] /= size;
+    }
+#endif
 
 	t = clock();	
 	// y[n] = sum(x[k], k=n-M..n+M)/(2M+1), n=0..N-1 
 	// zero padding for x[k] if k<0 or k>=N
-	for(n=0;n<N;n++){
-		y[n] = 0.0;
-		for(k=n-M;k<=n+M;++k){
-			if(k>=0 || k<N)
-				y[n] = y[n] + x[k];
-		}
-		y[n] = y[n]/(2*M+1);
-	}
+    //initialize y[0]
+    for(int i = 0; i <= M; ++i) y[0] += x[i];
+
+
+    // Use branchless method to reduce if-else
+    for(int i = 1, lower = - M, upper = 1 + M; i < N; ++i, ++lower, ++upper){
+        y[i] = (y[i - 1] - x[lower]) * (lower >= 0) + y[i-1]*(lower < 0);
+        y[i] += x[i + M] * (upper < N);
+    }
+
+    for(int i = 0; i < N; ++i) y[i] /= size;
+
 
 	// Output the results
 	t = clock() - t;
-	printf("y[0]=%f, y[1]=%f\n",y[0],y[1]);
-
+    for(int i = 0; i < N; ++i){
+#ifdef TEST
+        assert(y[i] - y_test[i] < 0.001); // 0.001 for floating number bias
+#endif
+        printf("%.3f\n", y[i]);
+    }
+	printf("\n");
 	// free the memory located by x, y
-
-
+    free(x);
+    free(y);
+
+#ifdef TEST
+    free(y_test);
+#endif
 
 	return 100;
 }

5_practice.c

@@ -2,37 +2,51 @@
 #include <stdlib.h>	// for memory allocation
 #include <time.h>	// for time calculation
 #include <math.h>	// for sine and cosine functions
+#include <string.h> // for memset
+
+
 int main() {
 	// Declare all the variables
-	int k, n, N;
+	int N;
 	double *x, *yr, *yi;
-	time_t t;
+	clock_t t;
 
 	// Input the number N
-
+    printf("Please input N = ");
+    scanf("%d", &N);    	
 
 	// Locate the memory for x, yr, yi;
+    x = (double*)malloc(sizeof(double)*N);
+    yr = (double*)malloc(sizeof(double)*N);
+    yi = (double*)malloc(sizeof(double)*N);
 
 	// Initial setting for x, for example, x[k] = k
-
-
-
-
+    for(int i = 0; i < N; ++i) x[i] = i;	
+
+    memset(yr, 0, sizeof(double)*N);
+    memset(yi, 0, sizeof(double)*N);
+
+    double pi_2_n_div_N, theta;
 	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 
-
-
-
-
-
+    for(int n = 0; n < N; ++n){
+        pi_2_n_div_N = 2*M_PI*n / (double)N;
+        for(int k = 0; k < N; ++k){
+            theta = pi_2_n_div_N * k;
+            yr[n] += cos(theta) * x[k];
+            yi[n] -= sin(theta) * x[k];
+        }
+    }
 
 
 	// output the results
 	t = clock() - t;
-	printf("%d ms for discrete Fourier Transform of %d elements\n", t, N);
+	printf("%.0f ms for discrete Fourier Transform of %d elements\n", 1000*(double)t/CLOCKS_PER_SEC, N);
 
 	// free the memory located by x, yr, yi
-
+    free(x);
+    free(yr);
+    free(yi);
 
 
 	return 100;

hw6/main.c

@@ -0,0 +1,142 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+
+
+struct __array_t{
+    unsigned int size;
+    unsigned int capacity;
+    double *re;
+    double *im;
+};
+typedef struct __array_t array_t;
+
+#define ARRAY_INIT (array_t){.size = 0, .capacity = 0, .re=NULL, .im=NULL}
+
+void push_data(array_t * ar, int data){
+    if(ar == NULL){
+        return; 
+    }
+
+    if(ar->size + 1 >= ar->capacity){
+        ar->re = realloc(ar->re, sizeof(double)*(ar->capacity + 1000));
+        ar->im = realloc(ar->im, sizeof(double)*(ar->capacity + 1000));
+        ar->capacity += 1000;
+    }
+
+    if(ar->re && ar->im){
+        ar->re[ar->size] = data;    
+        ++ar->size;
+    }
+}
+
+void prealloc_array(array_t *ar, unsigned int num){
+    ar->capacity = num;
+    ar->re = malloc(sizeof(double)*num);
+    ar->im = malloc(sizeof(double)*num);
+
+    if(!ar->re || !ar->im){
+        printf("Malloc failed in prealloc_array function\n");
+        exit(-1);
+    }
+}
+
+void free_array(array_t ar){
+    if(ar.re && ar.im){
+        free(ar.re);
+        free(ar.im);
+    }
+}
+
+void FFT2(double *xre, double *xim, double *yre, double *yim){
+    // y[0] = x[0] + x[1];
+    yre[0] = xre[0] + xre[1];
+    yim[0] = xim[0] + xim[1];
+
+    // y[1] = x[0] - x[1];
+    yre[1] = xre[0] - xre[1];
+    yim[1] = xim[0] - xim[1];
+}
+
+void FFTN(array_t x, array_t y, int N){
+
+    double cos_value, sin_value;
+
+    if(N & 2){
+        FFT2(x.re, x.im, y.re, y.im);
+    }else{
+        array_t xodd, xeven;
+        array_t yodd, yeven;
+        prealloc_array(&xodd, N>>1l);
+        prealloc_array(&xeven, N>>1l);
+
+        prealloc_array(&yodd, N>>1l);
+        prealloc_array(&yeven, N>>1l);
+
+
+        // place the data
+        for(int i = 0; i < N/2; ++i){
+            xodd.re[i] = x.re[(i<<1) | 1];
+            xodd.im[i] = x.im[(i<<1) | 1];
+
+            xeven.re[i] = x.re[i<<1];
+            xeven.im[i] = x.im[i<<1];
+        }
+
+        FFTN(xeven, yeven, N>>1l);
+        FFTN(xodd, yodd, N>>1l);
+
+        // collect the data
+        for(int i = 0; i <N/2; ++i){
+            cos_value = cos(2*M_PI*i/N);
+            sin_value = -sin(2*M_PI*i/N);
+
+            // for the value k < N/2
+            y.re[i] = yeven.re[i] + (yodd.re[i] *cos_value - yodd.im[i] * sin_value);
+            y.im[i] = yeven.im[i] + (yodd.re[i] *sin_value + yodd.im[i] * cos_value);
+            // for the value k >= N/2
+            y.re[N/2 + i] = yeven.re[i] - (yodd.re[i] *cos_value - yodd.im[i] * sin_value);
+            y.im[N/2 + i] = yeven.im[i] + (yodd.re[i] *sin_value + yodd.im[i] * cos_value);
+        }
+
+        free_array(xeven);
+        free_array(xodd);
+        free_array(yeven);
+        free_array(yodd);
+    }
+}
+
+
+int main(){
+    FILE * file = fopen("hw6.csv", "r");
+    if(file == NULL){
+        perror("Can't open file hw6.csv");
+        exit(-1);
+    }
+
+    array_t array = ARRAY_INIT;
+
+    int index, data;
+    while(fscanf(file,"%d,%d", &index, &data) != EOF){
+        push_data(&array, data);
+    }
+    fclose(file);
+
+    for(int i = 0; i < array.size; ++i){
+        printf("%d,%.2f\n", i, array.re[i]); 
+    }
+
+    free_array(array);
+
+
+    // truncate the data
+    array.size = 1l << (31l - __builtin_clz(array.size));
+    array_t array_res;
+    prealloc_array(&array_res, array.size);
+
+    FFTN(array, array_res, array.size);
+    array_res.size = array.size; 
+    for(int i = 0; i < array_res.size; ++i){
+        printf("%d,%.3f%s%.3f\n", i, array_res.re[i], (array_res.im[i]>0) ? "+": "", array_res.im[i]);
+    }
+}