diff --git a/tests/gensin1.c b/tests/gensin1.c new file mode 100644 index 0000000..81ca038 --- /dev/null +++ b/tests/gensin1.c @@ -0,0 +1,86 @@ +/* + * + * Copyright 2022 Oleg Borodin + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, + * MA 02110-1301, USA. + * + */ + + +#include +#include +#include +#include + + + +typedef struct { + double x0; + double x1; + double rc; +} lpf_t; + +void lpf_init(lpf_t *lpf, double freq) { + lpf->x0 = 0.0; + lpf->x1 = 0.0; + + double order = 1.5; + double n = 1 / sqrt(pow(2, 1.0 / order) - 1); + lpf->rc = 1 / (2 * n * M_PI * freq); +} + + +double lpf_next(lpf_t *lpf, double x, double dt) { + + double k = dt / (lpf->rc + dt); + + lpf->x1 = lpf->x1 + k * (x - lpf->x1); + lpf->x0 = lpf->x0 + k * (lpf->x1 - lpf->x0); + return lpf->x0; +} + +int main(int argc, char **argv) { + + double t = 0; + double s = 0; + double freq1 = 10; + double freq2 = 50; + double freq3 = 200; + double dt = 0.0001; + double l = 0.5; + + int fd = open("freq.dat", O_WRONLY | O_TRUNC | O_CREAT, 0666); + printf("start\n"); + + if (fd < 0) { return 1; }; + + lpf_t lpf; + lpf_init(&lpf, 30.0); + + + for (int i = 0; i < (l / dt); i++) { + t = i * dt; + double s1 = sin(2 * M_PI * freq1 * t); + double s2 = cos(2 * M_PI * freq2 * t); + double s3 = sin(2 * M_PI * freq3 * t); + s = (s1 + s2 + s3) / 3.0; + s = lpf_next(&lpf, s, dt); + dprintf(fd, "%e\t%e\n", t, s); + } + fsync(fd); + close(fd); + return 0; +} diff --git a/tests/gensin3.c b/tests/gensin3.c new file mode 100644 index 0000000..03c4a93 --- /dev/null +++ b/tests/gensin3.c @@ -0,0 +1,185 @@ +/* + * Copyright 2022 Oleg Borodin + */ + + +#include +#include +#include +#include +#include +#include + +typedef struct { + double x0; + double x1; + double rc; +} lpf_t; + +void lpf_init(lpf_t * lpf, double freq) { + lpf->x0 = 0.0; + lpf->x1 = 0.0; + + double order = 1.5; + double n = 1 / sqrt(pow(2, 1.0 / order) - 1); + + lpf->rc = 1 / (2 * n * M_PI * freq); +} + +double lpf_next(lpf_t * lpf, double x, double dt) { + + double k = dt / (lpf->rc + dt); + + lpf->x1 = lpf->x1 + k * (x - lpf->x1); + lpf->x0 = lpf->x0 + k * (lpf->x1 - lpf->x0); + return lpf->x0; +} + +void dft(double* x1, double* y1, int m, int dir) { + long i, k; + double arg; + double cosarg, sinarg; + + double x2[m]; + double y2[m]; + + for (i = 0; i < m; i++) { + x2[i] = 0; + y2[i] = 0; + arg = -dir * 2.0 * M_PI * (double)i / (double)m; + for (k = 0; k < m; k++) { + cosarg = cos(k * arg); + sinarg = sin(k * arg); + x2[i] += (x1[k] * cosarg - y1[k] * sinarg); + y2[i] += (x1[k] * sinarg + y1[k] * cosarg); + } + } + if (dir == 1) { + for (i = 0; i < m; i++) { + x1[i] = x2[i] / (double)m; + y1[i] = y2[i] / (double)m; + } + } else { + for (i = 0; i < m; i++) { + x1[i] = x2[i]; + y1[i] = y2[i]; + } + } +} + + +int main(int argc, char** argv) { + + double t = 0; + double s = 0; + double freq1 = 20; + double freq2 = 70; + double freq3 = 150; + double dt = 0.001; + + lpf_t lpf; + + lpf_init(&lpf, 20.0); + + int c = 1024 * 4; + double xa[c]; + double xb[c]; + double y[c]; + + for (int i = 0; i < c; i++) { + t = i * dt; + double s1 = sin(2 * M_PI * freq1 * t); + double s2 = cos(2 * M_PI * freq2 * t); + double s3 = sin(2 * M_PI * freq3 * t); + + s = (s1 + s2 + s3) / 3.0; + xa[i] = s; + xb[i] = lpf_next(&lpf, s, dt);; + } + + for (int i = 0; i < c; i++) { + y[i] = 0; + } + + double rate = 1.0 / dt; + + + FILE * fd; + + fd = popen ("gnuplot -persistent", "w"); + + fprintf(fd, "set terminal png size 1600,1200\n"); + fprintf(fd, "set output \"gensin3a.png\"\n"); + fprintf(fd, "plot '-' with lines\n"); + + dft(xa, y, c, -1); + for (int i = 0; i < c; i++) { + double freq = rate * (double)i / (double)c; + double amp = sqrt(xa[i]*xa[i] + y[i]*y[i]) ; + if (freq < rate / 2) { + fprintf(fd, "%e %e\n", freq, amp); + } + } + fflush(fd); + fprintf(fd, "e"); + fclose(fd); + + for (int i = 0; i < c; i++) { + y[i] = 0; + } + + + fd = popen ("gnuplot -persistent", "w"); + + fprintf(fd, "set terminal png size 1600,1200\n"); + fprintf(fd, "set output \"gensin3b.png\"\n"); + fprintf(fd, "plot '-' with lines\n"); + + dft(xb, y, c, -1); + for (int i = 0; i < c; i++) { + double freq = rate * (double)i / (double)c; + double amp = sqrt(xb[i]*xb[i] + y[i]*y[i]) ; + if (freq < rate / 2) { + fprintf(fd, "%e %e\n", freq, amp); + } + } + fflush(fd); + fprintf(fd, "e"); + fclose(fd); + + + +#if 0 + FILE* fd = fopen("freq1.dat", "w"); + + dft(xa, y, c, -1); + for (int i = 0; i < c; i++) { + double freq = rate * (double)i / (double)c; + double amp = sqrt(xa[i]*xa[i] + y[i]*y[i]) ; + if (freq < rate / 2) { + fprintf(fd, "%e %e\n", freq, amp); + } + } + + fclose(fd); + + for (int i = 0; i < c; i++) { + y[i] = 0; + } + + fd = fopen("freq2.dat", "w"); + + dft(xb, y, c, -1); + for (int i = 0; i < c; i++) { + double freq = rate * (double)i / (double)c; + double amp = sqrt(xb[i]*xb[i] + y[i]*y[i]); + if (freq < rate / 2) { + fprintf(fd, "%e %e\n", freq, amp); + } + } + + fclose(fd); +#endif + + return 0; +} diff --git a/tests/gensin3a.png b/tests/gensin3a.png new file mode 100644 index 0000000..6f41361 Binary files /dev/null and b/tests/gensin3a.png differ diff --git a/tests/gensin3b.png b/tests/gensin3b.png new file mode 100644 index 0000000..34da89a Binary files /dev/null and b/tests/gensin3b.png differ