diff --git a/.gitignore b/.gitignore index 321fb93..5d30fb2 100644 --- a/.gitignore +++ b/.gitignore @@ -1,8 +1,2 @@ fft_calc -sin01_fft.txt -sin02_fft.txt -sin03_fft.txt -sin04_fft.txt -sin800_fft.txt -sin2367_fft.txt - +*_fft* diff --git a/README.md b/README.md index 7870b88..9455291 100644 --- a/README.md +++ b/README.md @@ -19,8 +19,11 @@ except sin03 which is -1.1...1.1. * sin2367 - 4s sine wave 2367 Hz # Build -g++ -o fft_calc main.c fft.c +`g++ -o fft_calc main.cpp fft.c` # Usage -./fft_calc signal/sin01 signal/sin02 -./fft_calc signal/sin03 +`./fft_calc signal/sin01.txt signal/sin02.txt` \ +`./fft_calc signal/sin03.txt` + +To use the additional python script for plotting (requires python 3, matplotlib and numpy): \ +`python3 plots.py sin02` diff --git a/main.cpp b/main.cpp index ddac88c..a2bbd7e 100644 --- a/main.cpp +++ b/main.cpp @@ -73,12 +73,12 @@ int main(int argc, char *argv[]) { while (fgets(line, SIZE_OF_LINE, fp) != NULL) { - signal[i++] = atof(line); + signal[i++] = atof(line); } if (ferror(fp)) perror("Error: "); - printf("Read %u lines\n", i); + printf("Read %u lines\n", i); // entire signal read into memory fclose(fp); /***************************************************************************** @@ -210,6 +210,52 @@ int main(int argc, char *argv[]) { fclose(fp2); } + + strcpy(strstr(res_file_name, ".txt"), "_only_windows.txt\0"); + + fp2 = fopen(res_file_name, "w"); + + if (!fp2) + { + printf("Failed to open %s!\n", res_file_name); + } + else + { + printf("Writing results to file %s\n", res_file_name); + + // Print header row. + fprintf(fp2, "xaxis_fq\t"); + + m = calc_power(SIGNAL_BUF_SIZE); + for(k = 1; k <= m; k++) + { + if (k == m) + { + fprintf(fp2, "win%u\n", k); + } + else + { + fprintf(fp2, "win%u\t", k); + } + } + + // Print results to file. Full signal and each window result in separate column. + for (k = 0; k < MOVING_WINDOW_SIZE / 2; k++) + { + fprintf(fp2, "%lf\t", (fq_step_win * k)); + + for ( m = 0; m < num_windows - 1; m++) + { + fprintf(fp2, "%lf\t", fft[k + (MOVING_WINDOW_SIZE / 2) * m]); + } + + fprintf(fp2, "%lf\n", fft[k + (MOVING_WINDOW_SIZE / 2) * m]); + } + + fclose(fp2); + } + + } } diff --git a/plots.py b/plots.py new file mode 100644 index 0000000..1b03728 --- /dev/null +++ b/plots.py @@ -0,0 +1,43 @@ +import sys + +import matplotlib.pyplot as plt +import numpy as np + +# Signal filename can be passed as: +# python3 plots.py sin02 +if len(sys.argv) > 1: + file = sys.argv[1] +else: + file = "sin01" + +sampling_freq = 10000 # Hz +sampling_length = 6.5 #s + +s = np.loadtxt(f'./signals/{file}.txt') +t = np.linspace(0, sampling_length, len(s)) + +# Visualize original signal and fft of full signal +fft_res = np.loadtxt(f'./signals/{file}_fft.txt', skiprows=1, usecols=[0,1]) + +x_lim = -1 +x_axis_freqs = fft_res[0:x_lim, 0] +energy = fft_res[0:x_lim, 1] + +fig, (ax0, ax1, ax2) = plt.subplots(3, 1) +ax0.plot(t, s) +ax1.set_title("FFT on full signal", y=1.0, pad=-14) +ax1.plot(x_axis_freqs, energy) + +# Visualize fft of windows +fft_windows_res = np.loadtxt(f'./signals/{file}_fft_only_windows.txt', skiprows=1) + +# take mean of all windows for each frequency in x-axis +windows_mean = np.mean(fft_windows_res[0:-1, 1:-1], axis=1, dtype=np.float32) + +x_axis_freqs = fft_windows_res[0:x_lim, 0] +energy = np.transpose(windows_mean) + +ax2.set_title("FFT on windows (avg)", y=1.0, pad=-14) +ax2.plot(x_axis_freqs, energy) + +plt.show() \ No newline at end of file