From 47f43cf5a180544d09ecb42b8b6b255d619f4c29 Mon Sep 17 00:00:00 2001 From: Lino Mediavilla Date: Wed, 4 Aug 2021 13:01:20 +0200 Subject: [PATCH 1/3] fix readme --- README.md | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/README.md b/README.md index 7870b88..a9a39af 100644 --- a/README.md +++ b/README.md @@ -19,8 +19,8 @@ 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 From 13a2069c8ad4c9dcc7c9b0f1d14736e00b26118f Mon Sep 17 00:00:00 2001 From: Lino Mediavilla Date: Thu, 5 Aug 2021 18:43:09 +0200 Subject: [PATCH 2/3] Update readme & add plotting script --- README.md | 9 ++++++--- plots.py | 29 +++++++++++++++++++++++++++++ 2 files changed, 35 insertions(+), 3 deletions(-) create mode 100644 plots.py diff --git a/README.md b/README.md index a9a39af..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.cpp fft.c +`g++ -o fft_calc main.cpp fft.c` # Usage -./fft_calc signal/sin01.txt signal/sin02.txt \ -./fft_calc signal/sin03.txt +`./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/plots.py b/plots.py new file mode 100644 index 0000000..db2c87c --- /dev/null +++ b/plots.py @@ -0,0 +1,29 @@ +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)) + +fft_res = np.loadtxt(f'./signals/{file}_fft.txt', skiprows=1, usecols=[0,1]) +nyq = round(len(fft_res)/2) + +freq = fft_res[0:nyq, 0] +energy = fft_res[0:nyq, 1] + +fig, (ax0, ax1) = plt.subplots(2, 1) +ax0.plot(t, s) +ax1.plot(freq, energy) + +plt.show() \ No newline at end of file From 5fc49759bd531cfc4b5ed8e2d5dfbe0641f32756 Mon Sep 17 00:00:00 2001 From: Lino Mediavilla Date: Fri, 6 Aug 2021 12:55:39 +0200 Subject: [PATCH 3/3] Add fft-window visualization --- .gitignore | 8 +------- main.cpp | 50 ++++++++++++++++++++++++++++++++++++++++++++++++-- plots.py | 26 ++++++++++++++++++++------ 3 files changed, 69 insertions(+), 15 deletions(-) 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/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 index db2c87c..1b03728 100644 --- a/plots.py +++ b/plots.py @@ -3,7 +3,7 @@ import matplotlib.pyplot as plt import numpy as np -# signal filename can be passed as: +# Signal filename can be passed as: # python3 plots.py sin02 if len(sys.argv) > 1: file = sys.argv[1] @@ -16,14 +16,28 @@ 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]) -nyq = round(len(fft_res)/2) -freq = fft_res[0:nyq, 0] -energy = fft_res[0:nyq, 1] +x_lim = -1 +x_axis_freqs = fft_res[0:x_lim, 0] +energy = fft_res[0:x_lim, 1] -fig, (ax0, ax1) = plt.subplots(2, 1) +fig, (ax0, ax1, ax2) = plt.subplots(3, 1) ax0.plot(t, s) -ax1.plot(freq, energy) +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