diff --git a/5_moving_average.c b/5_moving_average.c index efffe2b..4b9800d 100644 --- a/5_moving_average.c +++ b/5_moving_average.c @@ -2,43 +2,87 @@ #include // for memory allocation #include // for time calculation #include // for sine and cosine functions -int main() { +#include + +// #define TEST + +#ifdef TEST +#include +#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= 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=0 || k= 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; } diff --git a/5_practice.c b/5_practice.c index 70d6651..8686c77 100644 --- a/5_practice.c +++ b/5_practice.c @@ -2,37 +2,51 @@ #include // for memory allocation #include // for time calculation #include // for sine and cosine functions +#include // 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; diff --git a/hw4 b/hw4 new file mode 160000 index 0000000..caaed10 --- /dev/null +++ b/hw4 @@ -0,0 +1 @@ +Subproject commit caaed102aa8dc22c7fbac05fa27d3597bd514fc0 diff --git a/hw5 .pdf b/hw5 .pdf new file mode 100644 index 0000000..f9354fa Binary files /dev/null and b/hw5 .pdf differ diff --git a/hw6/main.c b/hw6/main.c new file mode 100644 index 0000000..621e373 --- /dev/null +++ b/hw6/main.c @@ -0,0 +1,142 @@ +#include +#include +#include + + +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 + 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]); + } +}