Add ASK BPSK DQPSK(not implemented yet) FFT QPSK WAV

FFT/fft.c

@@ -0,0 +1,57 @@
+#include <math.h>
+#include <complex.h>
+#include <stdint.h>
+#include <stdlib.h>
+
+complex float* fft(complex float* p, uint32_t N) {
+    complex float w = cexpf(-2 * I * M_PI / N);
+
+    if (N == 1)
+        return p;
+
+    complex float* pe = (complex float*)malloc(sizeof(complex float) * (N / 2));
+    complex float* po = (complex float*)malloc(sizeof(complex float) * (N / 2));
+    for (int i = 0; i < N - 1; i += 2) {
+       pe[i / 2] = p[i];
+       po[i / 2] = p[i + 1];
+    }
+
+    complex float* xe = fft(pe, N / 2);
+    complex float* xo = fft(po, N / 2);
+
+    complex float* ps = (complex float*)malloc(sizeof(complex float) * N);
+
+    for (int i = 0; i < N / 2; i++) {
+        ps[i] = xe[i] + cpowf(w, (complex float)i) * xo[i];
+        ps[i + N / 2] = xe[i] - cpowf(w, (complex float)i) * xo[i];
+    }
+    
+    return ps;
+}
+
+
+complex float* ifft(complex float* p, uint32_t N) {
+    complex float w = ((complex float)1 / N) * cexpf(2 * I * M_PI / N);
+
+    if (N == 1)
+        return p;
+
+    complex float* pe = (complex float*)malloc(sizeof(complex float) * (N / 2));
+    complex float* po = (complex float*)malloc(sizeof(complex float) * (N / 2));
+    for (int i = 0; i < N - 1; i += 2) {
+       pe[i / 2] = p[i];
+       po[i / 2] = p[i + 1];
+    }
+
+    complex float* xe = fft(pe, N / 2);
+    complex float* xo = fft(po, N / 2);
+
+    complex float* ps = (complex float*)malloc(sizeof(complex float) * N);
+
+    for (int i = 0; i < N / 2; i++) {
+        ps[i] = xe[i] + cpowf(w, (complex float)i) * xo[i];
+        ps[i + N / 2] = xe[i] - cpowf(w, (complex float)i) * xo[i];
+    }
+    
+    return ps;
+}

FFT/fft.h

@@ -0,0 +1,8 @@
+#include <math.h>
+#include <complex.h>
+#include <stdint.h>
+#include <stdlib.h>
+
+complex float* fft(complex float* p, uint32_t N);
+
+complex float* ifft(complex float* p, uint32_t N);

FFT/fftwav.c

@@ -0,0 +1,87 @@
+#include <math.h>
+#include <stdio.h>
+#include <complex.h>
+#include <stdint.h>
+#include <stdlib.h>
+#include "../WAV/wav.h"
+
+complex float* fft(complex float* p, uint32_t N) {
+    complex float w = cexpf(-2 * I * M_PI / N);
+
+    if (N == 1)
+        return p;
+
+    complex float* pe = (complex float*)malloc(sizeof(complex float) * (N / 2));
+    complex float* po = (complex float*)malloc(sizeof(complex float) * (N / 2));
+    for (int i = 0; i < N - 1; i += 2) {
+       pe[i / 2] = p[i];
+       po[i / 2] = p[i + 1];
+    }
+
+    complex float* xe = fft(pe, N / 2);
+    complex float* xo = fft(po, N / 2);
+
+    complex float* ps = (complex float*)malloc(sizeof(complex float) * N);
+
+    for (int i = 0; i < N / 2; i++) {
+        ps[i] = xe[i] + cpowf(w, (complex float)i) * xo[i];
+        ps[i + N / 2] = xe[i] - cpowf(w, (complex float)i) * xo[i];
+    }
+    
+    return ps;
+}
+
+
+complex float* ifft(complex float* p, uint32_t N) {
+    complex float w = ((complex float)1 / N) * cexpf(2 * I * M_PI / N);
+
+    if (N == 1)
+        return p;
+
+    complex float* pe = (complex float*)malloc(sizeof(complex float) * (N / 2));
+    complex float* po = (complex float*)malloc(sizeof(complex float) * (N / 2));
+    for (int i = 0; i < N - 1; i += 2) {
+       pe[i / 2] = p[i];
+       po[i / 2] = p[i + 1];
+    }
+
+    complex float* xe = fft(pe, N / 2);
+    complex float* xo = fft(po, N / 2);
+
+    complex float* ps = (complex float*)malloc(sizeof(complex float) * N);
+
+    for (int i = 0; i < N / 2; i++) {
+        ps[i] = xe[i] + cpowf(w, (complex float)i) * xo[i];
+        ps[i + N / 2] = xe[i] - cpowf(w, (complex float)i) * xo[i];
+    }
+    
+    return ps;
+}
+
+int main() {
+    int wav_len;
+    double* wav_data = read_wav("s.wav", &wav_len);
+
+    int fsize = 1 << (int)floorf(log2f((float)wav_len));
+
+    complex float* cwav_data = (complex float*)malloc(sizeof(complex float) * fsize);
+
+    for (int i = 0; i < fsize; i++)
+        cwav_data[i] = (complex float)wav_data[i];
+
+    complex float* ft = fft(cwav_data, fsize);
+
+    FILE *gnuplot = popen("gnuplot -persistent", "w");
+    fprintf(gnuplot, "set title 'fft'\n");
+    fprintf(gnuplot, "set logscale x\n");
+    fprintf(gnuplot, "set logscale y\n");
+    fprintf(gnuplot, "plot '-' with linespoints title 'fft'\n");
+
+    for (int i = 0; i < fsize; i++) {
+        fprintf(gnuplot, "%f %f\n",((float)i / fsize) * (float)44100 ,cabsf(ft[i])); 
+    }
+    fprintf(gnuplot, "e\n");
+    pclose(gnuplot);
+
+    return 0;
+}