diff --git a/ASK/ask.c b/ASK/ask.c
new file mode 100644
index 0000000..f97d318
--- /dev/null
+++ b/ASK/ask.c
@@ -0,0 +1,109 @@
+#include <stdio.h>
+#include <math.h>
+#include <stdlib.h>
+#include "../WAV/wav.h"
+
+// Fréquence d'échantillonage
+#define FS 44100 
+// Temps d'échantillonage
+#define TS 0.05
+// Nombre d'échantillions
+#define NS (int)(FS*TS)
+// Amplitude du signal
+#define A 22767 
+// Fréquence de la porteuse
+#define FC 400 
+
+// Fonction de concaténation deux tableau
+double* concat_tab(double **t, int a, int b) {
+    double* tf = malloc(sizeof(double) * a * b);
+
+    int pos = 0;
+    for (int i = 0; i < b; i++) {
+        for (int j = 0; j < a; j++) {
+            tf[pos] = t[i][j];
+            pos++;
+        }
+    }
+
+    return tf;
+}
+
+// Fontion permetant la modulation d'un bit
+double* modulation_unique(int b) {
+    double* s = (double*)malloc(NS * sizeof(double));
+    
+    int c = (b == 0) ? A / 2 : A; 
+
+    for (int n = 0; n < NS; n++) {
+        s[n] = c * cos(2 * M_PI * FC * ((double)n / FS));
+    }
+    return s;
+}
+
+// Fonction permetant la modulation d'un suite de bits
+double* modulation(int* b, int nb_symbole) {
+    double** s = (double**)malloc(sizeof(double*) * nb_symbole);
+   
+    for (int k = 0; k < nb_symbole; k++) {
+        s[k] = modulation_unique(b[k]);
+    }
+
+    double *sf = concat_tab(s, NS, nb_symbole);
+
+    for (int i = 0; i < nb_symbole; i++) {
+        free(s[i]);
+    }
+    free(s);
+
+    return sf;
+}
+
+int* askdemodulation(double* r, int nb_symbole) {
+    int* sb = (int*)malloc(sizeof(int) * nb_symbole);
+
+    double thr = (A + A / 2.0) / 4.0;
+
+    for (int k = 0; k < nb_symbole; k++) {
+        double sum = 0;
+       for (int n = 0 ; n < NS; n++) {
+           sum += fabs(r[k * NS + n]);
+       }
+       double avg = sum / NS;
+       sb[k] = (avg > thr) ? 1 : 0;
+    }
+
+
+
+    return sb;
+}
+
+int main() {
+    int bits[] = {0, 1, 0, 1, 1, 1, 0, 1, 0, 0, 1, 0, 1, 1,0, 1, 0, 1, 1, 1, 0, 1, 0, 0, 1, 0, 1, 1,0, 1, 0, 1, 1, 1, 0, 1, 0, 0, 1, 0, 1, 1,0, 1, 0, 1, 1, 1, 0, 1, 0, 0, 1, 0, 1, 1};
+
+    int nb_symbole = (sizeof(bits) / sizeof(bits[0])) / 2;
+
+    for (int i = 0; i < 2 * nb_symbole; i++) {
+       printf("%d ", bits[i]);
+    }
+    printf("\n");
+
+    double* sf = modulation(bits, 2 * nb_symbole);
+
+    write_wav("signal.wav", sf, (double)NS * (double)nb_symbole);
+    
+    int wav_len;
+    double* wav_data = read_wav("signal.wav", &wav_len);
+
+    int* dbits = askdemodulation(wav_data, 2 * nb_symbole);
+
+    for (int i = 0; i < 2 * nb_symbole; i++) {
+       printf("%d ", dbits[i]);
+    }
+    printf("\n");
+
+    free(sf);
+    return 0;
+}
+
+
diff --git a/ASK/out b/ASK/out
new file mode 100755
index 0000000..3f88118
Binary files /dev/null and b/ASK/out differ
diff --git a/ASK/signal.wav b/ASK/signal.wav
new file mode 100644
index 0000000..bdec889
Binary files /dev/null and b/ASK/signal.wav differ
diff --git a/BPSK/bpsk.c b/BPSK/bpsk.c
new file mode 100644
index 0000000..3fcc3a5
--- /dev/null
+++ b/BPSK/bpsk.c
@@ -0,0 +1,199 @@
+#include <stdio.h>
+#include <math.h>
+#include <stdlib.h>
+
+// Fréquence d'échantillonage
+#define FS 44100 
+// Temps d'échantillonage
+#define TS 0.05
+// Nombre d'échantillions
+#define NS (FS*TS)
+// Amplitude du signal
+#define A 32767 
+// Fréquence de la porteuse
+#define FC 400 
+
+// Fonction de concaténation deux tableau
+double* concat_tab(double **t, int a, int b) {
+    double* tf = malloc(sizeof(double) * a * b);
+
+    int pos = 0;
+    for (int i = 0; i < b; i++) {
+        for (int j = 0; j < a; j++) {
+            tf[pos] = t[i][j];
+            pos++;
+        }
+    }
+
+    return tf;
+}
+
+// Fontion permetant la modulation d'un bit
+double* modulation_unique(int b) {
+    double* s = (double*)malloc(NS * sizeof(double));
+
+    // Phase de 0 pour un bit de 0 et phase de pi pour un bit de 1
+    int p = (b == 0) ? 0 : M_PI;
+    
+    for (int n = 0; n < NS; n++) {
+        s[n] = A * cos(2 * M_PI * FC * ((double)n / (double)FS) + p);
+    }
+    return s;
+}
+
+// Fonction permetant la modulation d'un suite de bits
+double* modulation(int* b, int bs) {
+    double** s = (double**)malloc(sizeof(double*) * bs);
+    
+    for (int k = 0; k < bs; k++) {
+        //s[k] = (double*)malloc(sizeof(double) * NS);
+        s[k] = modulation_unique(b[k]);
+    }
+
+    //double* sf = (double*)malloc(sizeof(double) * NS * bs);
+    double *sf = concat_tab(s, NS, bs);
+
+    for (int i = 0; i < bs; i++) {
+        free(s[i]);
+    }
+    free(s);
+
+    return sf;
+}
+
+// Demodulation bpsk en connaissant le nombre d'échantillions et la fréquence de la porteuse (et donc la fréquence d'échantillonage)
+int* demodulation(double* r, int rsize) {
+    int bs = rsize; // A MODIFIER QUAND LE TRUC DU MAIN MARCHERA !
+    int* b = (int*)malloc(sizeof(int) * bs);
+
+    int k = 0, i = 0; 
+    while (i < NS * bs) {
+        double tmp = 0;
+        for (int n = 0; n < NS; n++) {
+            tmp += r[k * (int)NS + n] * cos(2 * M_PI * FC * ((double)n / (double)FS));
+        }
+        // Developper le produit de cos et appliquer un filtre passe bas pour trouver un equivalent de tmp
+        b[k] = (tmp > 0.0) ? 0 : 1;
+        i += NS;
+        k++;
+    }
+    return b;
+}
+
+// Fonction pour écrire un fichier WAV 16 bits PCM
+void write_wav(const char* filename, double* data, int len) {
+    FILE* f = fopen(filename, "wb");
+    if (!f) {
+        return;
+    }
+
+    int32_t chunk_size = 36 + len * 2;
+    int16_t audio_format = 1; // PCM
+    int16_t num_channels = 1;
+    int32_t sample_rate = FS;
+    int32_t byte_rate = FS * num_channels * 2;
+    int16_t block_align = num_channels * 2;
+    int16_t bits_per_sample = 16;
+    int32_t subchunk2_size = len * 2;
+
+    // En-tête WAV
+    fwrite("RIFF", 1, 4, f);
+    fwrite(&chunk_size, 4, 1, f);
+    fwrite("WAVE", 1, 4, f);
+    fwrite("fmt ", 1, 4, f);
+
+    int32_t subchunk1_size = 16;
+    fwrite(&subchunk1_size, 4, 1, f);
+    fwrite(&audio_format, 2, 1, f);
+    fwrite(&num_channels, 2, 1, f);
+    fwrite(&sample_rate, 4, 1, f);
+    fwrite(&byte_rate, 4, 1, f);
+    fwrite(&block_align, 2, 1, f);
+    fwrite(&bits_per_sample, 2, 1, f);
+
+    fwrite("data", 1, 4, f);
+    fwrite(&subchunk2_size, 4, 1, f);
+
+    for (int i = 0; i < len; i++) {
+        int16_t sample = (int16_t)(data[i]);
+        fwrite(&sample, sizeof(int16_t), 1, f);
+    }
+
+    fclose(f);
+}
+
+
+double* read_wav(const char* filename, int* out_len) {
+    FILE* f = fopen(filename, "rb");
+    if (!f) {
+        perror("Erreur ouverture WAV");
+        return NULL;
+    }
+
+    fseek(f, 0, SEEK_END);
+    long filesize = ftell(f);
+    fseek(f, 0, SEEK_SET);
+
+    fseek(f, 44, SEEK_SET);
+
+    int len = (filesize - 44) / 2;
+    double* data = malloc(sizeof(double) * len);
+    int16_t sample;
+
+    for (int i = 0; i < len; i++) {
+        fread(&sample, sizeof(int16_t), 1, f);
+        data[i] = (double)sample;
+    }
+
+    fclose(f);
+
+    *out_len = len;
+    return data;
+}
+
+int main() {
+    int bits[] = {0, 1, 0, 1, 1, 1, 0, 1, 0, 0, 1, 0, 1, 1,0, 1, 0, 1, 1, 1, 0, 1, 0, 0, 1, 0, 1, 1,0, 1, 0, 1, 1, 1, 0, 1, 0, 0, 1, 0, 1, 1,0, 1, 0, 1, 1, 1, 0, 1, 0, 0, 1, 0, 1, 1};
+
+    int bs = sizeof(bits) / sizeof(bits[0]);
+
+    for (int i = 0; i < bs; i++) {
+       printf("%d ", bits[i]);
+    }
+    printf("\n");
+
+    double* sf = modulation(bits, bs);
+
+    write_wav("signal.wav", sf, NS * bs);
+    
+    /*
+    FILE *f = fopen("signal.csv", "w");
+
+    if (!f) {
+        perror("Erreur ouverture fichier");
+        free(sf);
+        return 1;
+    }
+
+    for (int i = 0; i < NS * bs; i++) {
+        fprintf(f, "%d,%f\n", i, sf[i]); 
+    }
+
+    fclose(f);
+    */
+
+    // NE MARCHE PAS RENVOIE 1 !
+    //int rsize = (sizeof(sf) / (int)sizeof(sf[0]));
+    
+    int wav_len;
+    double* wav_data = read_wav("signal3.wav", &wav_len);
+
+    int* dbits = demodulation(wav_data, bs);
+
+    for (int i = 0; i < bs; i++) {
+       printf("%d ", dbits[i]);
+    }
+    printf("\n");
+
+    free(sf);
+    return 0;
+}
diff --git a/BPSK/out b/BPSK/out
new file mode 100755
index 0000000..a549684
Binary files /dev/null and b/BPSK/out differ
diff --git a/BPSK/signal.wav b/BPSK/signal.wav
new file mode 100644
index 0000000..390ae0d
Binary files /dev/null and b/BPSK/signal.wav differ
diff --git a/BPSK/signal2.wav b/BPSK/signal2.wav
new file mode 100644
index 0000000..b588a88
Binary files /dev/null and b/BPSK/signal2.wav differ
diff --git a/BPSK/signal3.wav b/BPSK/signal3.wav
new file mode 100644
index 0000000..c0f1fe2
Binary files /dev/null and b/BPSK/signal3.wav differ
diff --git a/BPSK/signal3_fixed.wav b/BPSK/signal3_fixed.wav
new file mode 100644
index 0000000..d796611
Binary files /dev/null and b/BPSK/signal3_fixed.wav differ
diff --git a/FFT/fft.c b/FFT/fft.c
new file mode 100644
index 0000000..3851909
--- /dev/null
+++ b/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;
+}
diff --git a/FFT/fft.h b/FFT/fft.h
new file mode 100644
index 0000000..b9fff86
--- /dev/null
+++ b/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);
diff --git a/FFT/fftwav.c b/FFT/fftwav.c
new file mode 100644
index 0000000..9fe8db6
--- /dev/null
+++ b/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;
+}
diff --git a/FFT/out b/FFT/out
new file mode 100755
index 0000000..861e393
Binary files /dev/null and b/FFT/out differ
diff --git a/FFT/s.wav b/FFT/s.wav
new file mode 100644
index 0000000..bdec889
Binary files /dev/null and b/FFT/s.wav differ
diff --git a/FFT/st.wav b/FFT/st.wav
new file mode 100644
index 0000000..980b1b7
Binary files /dev/null and b/FFT/st.wav differ
diff --git a/FFT/st2.wav b/FFT/st2.wav
new file mode 100644
index 0000000..caf411f
Binary files /dev/null and b/FFT/st2.wav differ
diff --git a/QPSK/out b/QPSK/out
new file mode 100755
index 0000000..893a729
Binary files /dev/null and b/QPSK/out differ
diff --git a/QPSK/qpsk.c b/QPSK/qpsk.c
new file mode 100644
index 0000000..1a3057c
--- /dev/null
+++ b/QPSK/qpsk.c
@@ -0,0 +1,108 @@
+#include <stdio.h>
+#include <math.h>
+#include <stdlib.h>
+#include "../WAV/wav.h"
+
+// Fréquence d'échantillonage
+#define FS 44100 
+// Temps d'échantillonage
+#define TS 0.05
+// Nombre d'échantillions
+#define NS (int)(FS*TS)
+// Amplitude du signal
+#define A 22767 
+// Fréquence de la porteuse
+#define FC 400 
+
+// Fonction de concaténation deux tableau
+double* concat_tab(double **t, int a, int b) {
+    double* tf = malloc(sizeof(double) * a * b);
+
+    int pos = 0;
+    for (int i = 0; i < b; i++) {
+        for (int j = 0; j < a; j++) {
+            tf[pos] = t[i][j];
+            pos++;
+        }
+    }
+
+    return tf;
+}
+
+// Fontion permetant la modulation d'une paire de bits
+double* modulation_unique(int b0, int b1) {
+    double* s = (double*)malloc(NS * sizeof(double));
+
+    //double p = (M_PI / 2) * (2 * b0 + b1);
+    int eps0 = 2 * b0 - 1; // Negatif si b = 0 et positif si b = 1
+    int eps1 = 2 * b1 - 1; // Negatif si b = 0 et positif si b = 1
+    
+    for (int n = 0; n < NS; n++) {
+        s[n] = A * eps0 * cos(2 * M_PI * FC * ((double)n / FS)) + A * eps1 * sin(2 * M_PI * FC * ((double)n / FS));
+    }
+    return s;
+}
+
+// Fonction permetant la modulation d'un suite de bits
+double* modulation(int* b, int nb_symbole) {
+    double** s = (double**)malloc(sizeof(double*) * nb_symbole );
+   
+    for (int k = 0; k < 2 * nb_symbole - 1; k += 2) {
+        s[k / 2] = modulation_unique(b[k], b[k + 1]);
+    }
+
+    double *sf = concat_tab(s, NS, nb_symbole);
+
+    for (int i = 0; i < nb_symbole; i++) {
+        free(s[i]);
+    }
+    free(s);
+
+    return sf;
+}
+
+int* qpskdemodulation(double* r, int nb_symbole) {
+    int* sb = (int*)malloc(sizeof(int) * 2 * nb_symbole);
+
+    for (int k = 0; k < nb_symbole; k++) {
+       double I = 0;
+       double Q = 0;
+
+       for (int n = 0 ; n < NS; n++) {
+           I += r[k * (int)NS + n] * A * cos(2 * M_PI * FC * ((double)n / (double)FS));
+           Q += r[k * (int)NS + n] * A * sin(2 * M_PI * FC * ((double)n / (double)FS));
+       }
+       sb[2 * k] = (I > 0) ? 1 : 0;
+       sb[2 * k + 1] = (Q > 0) ? 1 : 0;
+    }
+    return sb;
+}
+
+int main() {
+    int bits[] = {0, 1, 0, 1, 1, 1, 0, 1, 0, 0, 1, 0, 1, 1,0, 1, 0, 1, 1, 1, 0, 1, 0, 0, 1, 0, 1, 1,0, 1, 0, 1, 1, 1, 0, 1, 0, 0, 1, 0, 1, 1,0, 1, 0, 1, 1, 1, 0, 1, 0, 0, 1, 0, 1, 1};
+
+    int nb_symbole = (sizeof(bits) / sizeof(bits[0])) / 2;
+
+    for (int i = 0; i < 2 * nb_symbole; i++) {
+       printf("%d ", bits[i]);
+    }
+    printf("\n");
+
+    double* sf = modulation(bits, nb_symbole);
+
+    write_wav("signal.wav", sf, (double)NS * (double)nb_symbole);
+    
+    int wav_len;
+    double* wav_data = read_wav("signal.wav", &wav_len);
+
+    int* dbits = qpskdemodulation(wav_data, nb_symbole);
+
+    for (int i = 0; i < 2 * nb_symbole; i++) {
+       printf("%d ", dbits[i]);
+    }
+    printf("\n");
+
+    free(sf);
+    return 0;
+}
+
diff --git a/QPSK/signal.wav b/QPSK/signal.wav
new file mode 100644
index 0000000..caf411f
Binary files /dev/null and b/QPSK/signal.wav differ
diff --git a/README.md b/README.md
index e69de29..a3ef240 100644
--- a/README.md
+++ b/README.md
@@ -0,0 +1 @@
+DSP related implementations.
diff --git a/WAV/wav.c b/WAV/wav.c
new file mode 100644
index 0000000..199dc9e
--- /dev/null
+++ b/WAV/wav.c
@@ -0,0 +1,76 @@
+#include <stdio.h>
+#include <math.h>
+#include <stdlib.h>
+
+#define FS 44100
+
+// Fonction pour écrire un fichier WAV 16 bits PCM
+void write_wav(const char* filename, double* data, int len) {
+    FILE* f = fopen(filename, "wb");
+    if (!f) {
+        return;
+    }
+
+    int32_t chunk_size = 36 + len * 2;
+    int16_t audio_format = 1; // PCM
+    int16_t num_channels = 1;
+    int32_t sample_rate = FS;
+    int32_t byte_rate = FS * num_channels * 2;
+    int16_t block_align = num_channels * 2;
+    int16_t bits_per_sample = 16;
+    int32_t subchunk2_size = len * 2;
+
+    // En-tête WAV
+    fwrite("RIFF", 1, 4, f);
+    fwrite(&chunk_size, 4, 1, f);
+    fwrite("WAVE", 1, 4, f);
+    fwrite("fmt ", 1, 4, f);
+
+    int32_t subchunk1_size = 16;
+    fwrite(&subchunk1_size, 4, 1, f);
+    fwrite(&audio_format, 2, 1, f);
+    fwrite(&num_channels, 2, 1, f);
+    fwrite(&sample_rate, 4, 1, f);
+    fwrite(&byte_rate, 4, 1, f);
+    fwrite(&block_align, 2, 1, f);
+    fwrite(&bits_per_sample, 2, 1, f);
+
+    fwrite("data", 1, 4, f);
+    fwrite(&subchunk2_size, 4, 1, f);
+
+    for (int i = 0; i < len; i++) {
+        int16_t sample = (int16_t)(data[i]);
+        fwrite(&sample, sizeof(int16_t), 1, f);
+    }
+
+    fclose(f);
+}
+
+
+double* read_wav(const char* filename, int* out_len) {
+    FILE* f = fopen(filename, "rb");
+    if (!f) {
+        perror("Fail to open .wav");
+        return NULL;
+    }
+
+    fseek(f, 0, SEEK_END);
+    long filesize = ftell(f);
+    fseek(f, 0, SEEK_SET);
+
+    fseek(f, 44, SEEK_SET);
+
+    int len = (filesize - 44) / 2;
+    double* data = malloc(sizeof(double) * len);
+    int16_t sample;
+
+    for (int i = 0; i < len; i++) {
+        fread(&sample, sizeof(int16_t), 1, f);
+        data[i] = (double)sample;
+    }
+
+    fclose(f);
+
+    *out_len = len;
+    return data;
+}
diff --git a/WAV/wav.h b/WAV/wav.h
new file mode 100644
index 0000000..125cdda
--- /dev/null
+++ b/WAV/wav.h
@@ -0,0 +1,3 @@
+void write_wav(const char* filename, double* data, int len);
+
+double* read_wav(const char* filename, int* out_len);