小波變換C代碼

#include <stdio.h>
#include <stdlib.h>
#define LENGTH 512//信號長度
/******************************************************************
*　一維卷積函數
*
*　說明: 循環卷積,卷積結果的長度與輸入信號的長度相同
*
*　輸入參數: data[],輸入信號; core[],卷積核; cov[],卷積結果;
*             n,輸入信號長度; m,卷積核長度.
*
*　李承宇, lichengyu2345@126.com
*
*  2010-08-18  
******************************************************************/
void Covlution(double data[], double core[], double cov[], int n, int m)
{
    int i = 0;
    int j = 0;
    int k = 0;

    //將cov[]清零
    for(i = 0; i < n; i++)
    {
        cov[i] = 0;
    }

    //前m/2+1行
    i = 0;
    for(j = 0; j < m/2; j++, i++)
    {
        for(k = m/2-j; k < m; k++ )
        {
            cov[i] += data[k-(m/2-j)] * core[k];//k針對core[k]
        }

        for(k = n-m/2+j; k < n; k++ )
        {
            cov[i] += data[k] * core[k-(n-m/2+j)];//k針對data[k]
        }
    }

    //中間的n-m行
    for( i = m/2; i <= (n-m)+m/2; i++)
    {
        for( j = 0; j < m; j++)
        {
            cov[i] += data[i-m/2+j] * core[j];
        }
    }

    //最後m/2-1行
    i = (n - m) + m/2 + 1;
    for(j = 1; j < m/2; j++, i++)
    {
        for(k = 0; k < j; k++)
        {
            cov[i] += data[k] * core[m-j-k];//k針對data[k]
        }

        for(k = 0; k < m-j; k++)
        {
            cov[i] += core[k] * data[n-(m-j)+k];//k針對core[k]
        }
    }

}

/******************************************************************
*　一維小波變換函數
*
*　說明: 一維小波變換,只變換一次
*
*　輸入參數: input[],輸入信號; output[],小波變換結果，包括尺度係數和
*　小波係數兩部分; temp[],存放中間結果;h[],Daubechies小波基低通濾波器係數;
*　g[],Daubechies小波基高通濾波器係數;n,輸入信號長度; m,Daubechies小波基緊支集長度.
*
*　李承宇, lichengyu2345@126.com
*
*  2010-08-19  
******************************************************************/
void DWT1D(double input[], double output[], double temp[], double h[], 
           double g[], int n, int m)
{
//    double temp[LENGTH] = {0};//?????????????
    
    int i = 0;
/*
    //尺度係數和小波係數放在一塊兒
    Covlution(input, h, temp, n, m);

    for(i = 0; i < n; i += 2)
    {
        output[i] = temp[i];
    }

    Covlution(input, g, temp, n, m);

    for(i = 1; i < n; i += 2)
    {
        output[i] = temp[i];
    }
*/

    //尺度係數和小波係數分開
    Covlution(input, h, temp, n, m);

    for(i = 0; i < n; i += 2)
    {
        output[i/2] = temp[i];//尺度係數
    }

    Covlution(input, g, temp, n, m);

    for(i = 1; i < n; i += 2)
    {
        output[n/2+i/2] = temp[i];//小波係數
    }

}

void main()
{

    double data[LENGTH];//輸入信號
    double temp[LENGTH];//中間結果
    double data_output[LENGTH];//一維小波變換後的結果
    int n = 0;//輸入信號長度
    int m = 6;//Daubechies正交小波基長度
    int i = 0; 
    char s[32];//從txt文件中讀取一行數據

    static double h[] = {.332670552950, .806891509311, .459877502118, -.135011020010, 
                    -.085441273882, .035226291882};
    static double g[] = {.035226291882, .085441273882, -.135011020010, -.459877502118,
                    .806891509311, -.332670552950};
    //讀取輸入信號
    FILE *fp;
    fp=fopen("data.txt","r");
    if(fp==NULL) //若是讀取失敗
    {
        printf("錯誤！找不到要讀取的文件/"data.txt/"/n");
        exit(1);//停止程序
    }

    while( fgets(s, 32, fp) != NULL )//讀取長度n要設置得長一點，要保證讀到回車符，這樣指針纔會定位到下一行？回車符返回的是零值？是，非數字字符通過atoi變換都應該返回零值
    {
    //    fscanf(fp,"%d", &data[count]);//必定要有"&"啊！！！最後讀了個回車符！適應能力不如atoi啊
        data[n] = atof(s);
        n++;
    }

    //一維小波變換
    DWT1D(data, data_output, temp, h, g, n, m);

    //一維小波變換後的結果寫入txt文件
    fp=fopen("data_output.txt","w");

    //打印一維小波變換後的結果
    for(i = 0; i < n; i++)
    {
        printf("%f/n", data_output[i]);
        fprintf(fp,"%f/n", data_output[i]);
         }
    //關閉文件
    fclose(fp);
}