Tuesday, April 1, 2014

FSK 實驗筆記

1. 從 Wiki 下載範例檔案 http://en.wikipedia.org/wiki/File:AFSK_1200_baud.ogg

轉成 WAV 檔, 波形大概是這樣



位元被編碼成 1200 Hz (mark, 代表 1) 與 2200 Hz (space, 代表 0) 的訊號

2. 從小抄貼一些計算過來

令編碼後的訊號為 x(n) = A * sin(2pi * F * n * T)

x(n - k) 為 x(n) 延遲 k 個樣本的訊號

v(n) = x(n) * x(n - k)

= A^2 * sin(2pi * F * n * T) * sin(2pi * F * (n - k) * T)

= A^2 * [cos(2pi * F * k * T) - cos(4pi * F * n * T - 2pi * F * k * T)] / 2

後項因為頻率較高, 可以用低通濾波去掉

而前項已跟 n 無關, 依現在訊號代表的位元呈現兩種不同的常數

3. 計算濾波器參數

說是計算, 但是現在應該沒有人吃飽太閒用手算了

用 Scilab 生一個 5 階的, 截止頻率約 1200 Hz 的 Butterworth 低通濾波器

(小抄說 3 階就夠, 我怎麼轉出來有點抖)


0.027211 ~= 1200 / 44100 (取樣率)

4. 轉換 WAV 檔

上面才剛說, 這邊突然又想不開用 C 了

#include <stdio.h>
#include <string.h>

// assume the WAV header length is 44 bytes
#define HEADER_LEN       44
#define BUF_LEN          4096
#define DELAY            4

#define ARRAY_LEN(a)     (sizeof(a) / sizeof(a[0]))
#define ARRAY_SHIFT(a)   memmove(a + 1, a, (ARRAY_LEN(a) - 1)*sizeof(a[0]))

// assume the WAV file format is 44100 Hz/16 bit/Mono
// the cutoff frequency is about 1200 Hz
#define LP_ORDER         5
#define LP_A_COEF        1, -4.4469186, 7.9375491, -7.1066403, 3.1907047, -0.5745827
#define LP_B_COEF        0.0000035, 0.0000175, 0.0000350, 0.0000350, 0.0000175, 0.0000035

double lp(int v[])
{
    static const double a[] = {LP_A_COEF};
    static const double b[] = {LP_B_COEF};
    static double y[LP_ORDER + 1] = {0};
    double output;
    int i;

    y[0] = b[0]*v[0];

    for (i = 1; i <= LP_ORDER; i++)
        y[0] += b[i]*v[i] - a[i]*y[i];

    output = y[0];
    ARRAY_SHIFT(y);
    return output;
}

int main(int argc, char *argv[])
{
    char header[HEADER_LEN];
    short in[BUF_LEN];
    short out[BUF_LEN];
    int buf_usage;
    int x[DELAY + 1] = {0};
    int v[LP_ORDER + 1] = {0};
    FILE *file_in;
    FILE *file_out;

    if (argc != 3)
        return 1;

    file_in = fopen(argv[1], "rb");
    file_out = fopen(argv[2], "wb");

    if (!file_in || !file_out)
        return 2;

    fread(header, 1, HEADER_LEN, file_in);
    fwrite(header, 1, HEADER_LEN, file_out);

    while ((buf_usage = fread(in, sizeof(short), BUF_LEN, file_in)) != 0) {
        int i;

        for (i = 0; i < buf_usage; i++) {
            x[0] = in[i];
            v[0] = (x[0]*x[DELAY]) >> 15;
            out[i] = (short)lp(v);
            ARRAY_SHIFT(x);
            ARRAY_SHIFT(v);
        }

        fwrite(out, sizeof(short), buf_usage, file_out);
    }

    fclose(file_in);
    fclose(file_out);
    return 0;
}

5. 輸出結果


根據振幅可得知一段區間內所代表的位元

掰不下去了, 邊打邊看中天邊笑

記者都不會不好意思喔

小抄: Teaching DSP through the Practical Case Study of an FSK Modem