/**
 * @brief 小波变换之分解
 * @param sourceData 源数据
 * @param dataLen 源数据长度
 * @param db 过滤器类型
 * @param cA 分解后的近似部分序列-低频部分
 * @param cD 分解后的细节部分序列-高频部分
 * @return
 *         正常则返回分解后序列的数据长度，错误则返回-1
 */ 
int Wavelet::Dwt(double *sourceData, int dataLen, Filter db, double *cA, double *cD) 
{ 
    if(dataLen < 2) 
        return -1; 
    if((NULL == sourceData)||(NULL == cA)||(NULL == cD)) 
        return -1; 
 
    m_db = db; 
    int filterLen = m_db.length; 
    int n,k,p; 
    int decLen = (dataLen+filterLen-1)/2; 
    double tmp=0; 
    cout<<"decLen="<<decLen<<endl; 
 
    for(n=0; n<decLen; n++) 
    { 
        cA[n] = 0; 
        cD[n] = 0; 
        for(k=0; k<filterLen; k++) 
        { 
            p = 2*n-k;  // 感谢网友onetrain指正，此处由p = 2*n-k+1改为p = 2*n-k，解决小波重构后边界异常问题--2013.3.29 by hmm 
 
            // 信号边沿对称延拓 
            if((p<0)&&(p>=-filterLen+1)) 
                tmp = sourceData[-p-1]; 
            else if((p>dataLen-1)&&(p<=dataLen+filterLen-2)) 
                tmp = sourceData[2*dataLen-p-1]; 
            else if((p>=0)&&(p<dataLen-1)) 
                tmp = sourceData[p]; 
            else 
                tmp = 0; 
 
            // 分解后的近似部分序列-低频部分 
            cA[n] += m_db.lowFilterDec[k]*tmp; 
 
            // 分解后的细节部分序列-高频部分 
            cD[n] += m_db.highFilterDec[k]*tmp; 
        } 
 
    } 
 
    return decLen; 
} 
2.5      小波重构
[cpp] view plaincopyprint?
/**
 * @brief 小波变换之重构
 * @param cA 分解后的近似部分序列-低频部分
 * @param cD 分解后的细节部分序列-高频部分
 * @param cALength 输入数据长度
 * @param db 过滤器类型
 * @param recData 重构后输出的数据
 */ 
void  Wavelet::Idwt(double *cA, double *cD,  int cALength, Filter db, double *recData) 
{ 
    if((NULL == cA)||(NULL == cD)||(NULL == recData)) 
        return; 
 
    m_db = db; 
    int filterLen = m_db.length; 
 
    int n,k,p; 
    int recLen = 2*cALength-filterLen+1; 
    cout<<"recLen="<<recLen<<endl; 
 
    for(n=0; n<recLen; n++) 
    { 
        recData[n] = 0; 
        for(k=0; k<cALength; k++) 
        { 
            p = n-2*k+filterLen-1; 
 
            // 信号重构 
            if((p>=0)&&(p<filterLen)) 
            { 
                recData[n] += m_db.lowFilterRec[p]*cA[k] + m_db.highFilterRec[p]*cD[k]; 
                //cout<<"recData["<<n<<"]="<<recData[n]<<endl; 
            } 
 
        } 
    } 
} 
 
 
 
 
 
//***************************************************************************************************
  1 #define LENGTH 512
  2 #define LEVEL 4
  3 #define L_core 6
  4
  5 static void Covlution(double data[], double core[], double cov[], int LEN)
  6 {
  7         double temp[LENGTH + L_core - 1] = {0};
  8         int i = 0;
  9         int j = 0;
 10
 11     for(i = 0; i < LEN; i++)
 12         {
 13                 for(j = 0; j  < L_core; j++)
 14                 {
 15             temp[i + j] += data[i] * core[j];
 16         }
 17     }
 18
 19         for(i = 0; i < LEN; i++)
 20         {
 21                 if(i < L_core - 1)
 22                         cov[i] = temp[i] + temp[LEN + i];
 23                 else
 24                         cov[i] = temp[i];
 25         }
 26
 27 }
 28
 29 static void Covlution2(double data[], double core[], double cov[], int LEN)
 30 {
 31         double temp[LENGTH + L_core - 1] = {0};
 32         int i = 0;
 33         int j = 0;
 34
 35     for(i = 0; i < LEN; i++)
 36         {
 37                 for(j = 0; j  < L_core; j++)
 38                 {
 39             temp[i + j] += data[i] * core[j];
 40         }
 41     }
 42
 43         for(i = 0; i < LEN; i++)
 44         {
 45                 if(i < L_core - 1)
 46                         cov[i + LEN - L_core + 1] = temp[i] + temp[LEN + i];
 47                 else
 48                         cov[i - L_core + 1] = temp[i];
 49         }
 50
 51 }
 52
 53 static void DWT1D(double input[], double output[], double LF[], double HF[], int l)
 54 {
 55     int i = 0;
 56         double temp[LENGTH] = {0};
 57         int LEN = LENGTH / pow(2, l - 1);
 58
 59     Covlution(input, LF, temp, LEN);
 60     for(i = 1; i < LEN; i += 2)
 61     {
 62         output[i/2] = temp[i];
 63     }
 64
 65     Covlution(input, HF, temp, LEN);
 66     for(i = 1; i < LEN; i += 2)
 67     {
 68         output[LEN/2 + i/2] = temp[i];
 69     }
 70 }
 71
 72 static void DWT(double input[], double output[], double LF[], double HF[], int len[])
 73 {
 74         int i;
 75         int j;
 76
 77         len[0] = len[1] = LENGTH / pow(2, LEVEL);
 78         for(i = 2; i <= LEVEL; i++) len[i] = len[i - 1] * 2;
 79
 80         DWT1D(input, output, LF, HF, 1);
 81         for(i = 2; i <= LEVEL; i++)
 82     {
 83                 for(j = 0; j < len[LEVEL + 2 - i]; j++) input[j] = output[j];
 84         DWT1D(input, output, LF, HF, i);
 85     }
 86 }
 87
 88 static void IDWT1D(double input[], double output[], double LF[], double HF[], int l, int flag)
 89 {
 90     int i = 0;
 91         double temp[LENGTH] = {0};
 92         int LEN = l * 2;
 93
 94         if(flag) Covlution2(input, HF, temp, LEN);
 95         else Covlution2(input, LF, temp, LEN);
 96
 97         for(i = 0; i < LEN; i++)
 98     {
 99         output[i] = temp[i];
100     }
101 }
102
103 static void IDWT(double input[], double output[], double LF[], double HF[], int len[], int level)
104 {
105         int i;
106         int j;
107         for(j = 0; j < len[LEVEL + 1 - level]; j++)
108         {
109                 output[2 * j] = 0;
110                 output[2 * j + 1] = input[j];
111         }
112         for(j = 0; j < 2 * len[LEVEL + 1 - level]; j++)
113         {
114                 input[j] = output[j];
115         }
116     IDWT1D(input, output, LF, HF, len[LEVEL + 1 - level], 1);
117
118         for(i = level - 1; i > 0; i--)
119     {
120                 for(j = 0; j < len[LEVEL + 1 - i]; j++)
121                 {
122                         input[2 * j] = 0;
123                         input[2 * j + 1] = output[j];
124                 }
125         IDWT1D(input, output, LF, HF, len[LEVEL + 1 - i], 0);
126     }
127 }