Haar小波变换基本原理

另外参见俄罗斯写的http://www.codeproject.com/Articles/22243/Real-Time-Object-Tracker-in-C Haar小波在图像处理和数字水印等方面应用较多，这里简单的介绍一下哈尔小波的基本原理以及其实现情况。  一、Haar小波的基本原理         数学理论方面的东西我也不是很熟悉，这边主要用简单的例子来介绍下Haar小波的使用情况。          例如：有a=[8,7,6,9]四个数，并使用b[4]数组来保存结果.         则一级Haar小波变换的结果为:         b[0]=(a[0]+a[1])/2,                       b[2]=(a[0]-a[1])/2         b[1]=(a[2]+a[3])/2,                       b[3]=(a[2]-a[3])/2        即依次从数组中取两个数字，计算它们的和以及差，并将和一半和差的一半依次保存在数组的前半部分和后半部分。        例如：有a[8]，要进行一维Haar小波变换，结果保存在b[8]中         则一级Haar小波变换的结果为:         b[0]=(a[0]+a[1])/2,                        b[4]=(a[0]-a[1])/2         b[1]=(a[2]+a[3])/2,                        b[5]=(a[2]-a[3])/2         b[2]=(a[4]+a[5])/2,                        b[6]=(a[4-a[5]])/2         b[3]=(a[6]+a[7])/2,                        b[7]=(a[6]-a[7])/2                 如果需要进行二级Haar小波变换的时候，只需要对b[0]-b[3]进行Haar小波变换.         对于二维的矩阵来讲，每一级Haar小波变换需要先后进行水平方向和竖直方向上的两次一维小波变换,行和列的先后次序对结果不影响。 二、Haar小波的实现 使用opencv来读取图片及像素，对图像的第一个8*8的矩阵做了一级小波变换 #include <cv.h> #include <highgui.h> #include <iostream> using namespace std; int main(){ IplImage* srcImg; double  imgData[8][8]; int i,j;     srcImg=cvLoadImage(“lena.bmp”,0); cout<<”原8*8数据”<<endl; for( i=0;i<8;i++) { for( j=0;j<8;j++) { imgData[i][j]=cvGetReal2D(srcImg,i+256,j+16); cout<<imgData[i][j]<<” “; } cout<<endl; }     double tempData[8];     //行小波分解    for( i=0;i<8;i++)    {    for( j=0;j<4;j++)    { double temp1=imgData[i][2*j]; double temp2=imgData[i][2*j+1]; tempData[j]=(temp1+temp2)/2; tempData[j+4]=(temp1-temp2)/2;    } for( j=0;j<8;j++)  imgData[i][j]=tempData[j];     } //列小波分解 for( i=0;i<8;i++) { for( j=0;j<4;j++)    {    double temp1=imgData[2*j][i];    double temp2=imgData[2*j+1][i];    tempData[j]=(temp1+temp2)/2;        tempData[j+4]=(temp1-temp2)/2;    } for( j=0;j<8;j++)  imgData[j][i]=tempData[j]; } cout<<”1级小波分解数据”<<endl; for( i=0;i<8;i++) { for( j=0;j<8;j++) { cout<<imgData[i][j]<<” “; } cout<<endl; } //列小波逆分解 for( i=0;i<8;i++) { for( j=0;j<4;j++) { double temp1=imgData[j][i]; double temp2=imgData[j+4][i]; tempData[2*j]=temp1+temp2; tempData[2*j+1]=temp1-temp2; } for( j=0;j<8;j++) { imgData[j][i]=tempData[j]; } } //行小波逆分解 for( i=0;i<8;i++) { for( j=0;j<4;j++) { double temp1=imgData[i][j]; double temp2=imgData[i][j+4]; tempData[2*j]=temp1+temp2; tempData[2*j+1]=temp1-temp2; } for( j=0;j<2*4;j++) { imgData[i][j]=tempData[j]; } } cout<<”1级小波逆分解数据”<<endl; for( i=0;i<8;i++) { for( j=0;j<8;j++) { cout<<imgData[i][j]<<” “; } cout<<endl; }      return 0; } =========================================================================== /// 小波变换 Mat WDT( const Mat &_src, const string _wname, const int _level )const { int reValue = THID_ERR_NONE; Mat src = Mat_<float>(_src); Mat dst = Mat::zeros( src.rows, src.cols, src.type() ); int N = src.rows; int D = src.cols;   /// 高通低通滤波器 Mat lowFilter; Mat highFilter; wavelet( _wname, lowFilter, highFilter );   /// 小波变换 int t=1; int row = N; int col = D;   while( t<=_level ) { ///先进行行小波变换 for( int i=0; i<row; i++ ) { /// 取出src中要处理的数据的一行 Mat oneRow = Mat::zeros( 1,col, src.type() ); for ( int j=0; j<col; j++ ) { oneRow.at<float>(0,j) = src.at<float>(i,j); } oneRow = waveletDecompose( oneRow, lowFilter, highFilter ); /// 将src这一行置为oneRow中的数据 for ( int j=0; j<col; j++ ) { dst.at<float>(i,j) = oneRow.at<float>(0,j); } }   #if 0 //normalize( dst, dst, 0, 255, NORM_MINMAX ); IplImage dstImg1 = IplImage(dst); cvSaveImage( "dst.jpg", &dstImg1 ); #endif /// 小波列变换 for ( int j=0; j<col; j++ ) { /// 取出src数据的一行输入 Mat oneCol = Mat::zeros( row, 1, src.type() ); for ( int i=0; i<row; i++ ) { oneCol.at<float>(i,0) = dst.at<float>(i,j); } oneCol = ( waveletDecompose( oneCol.t(), lowFilter, highFilter ) ).t(); for ( int i=0; i<row; i++ ) { dst.at<float>(i,j) = oneCol.at<float>(i,0); } }   #if 0 //normalize( dst, dst, 0, 255, NORM_MINMAX ); IplImage dstImg2 = IplImage(dst); cvSaveImage( "dst.jpg", &dstImg2 ); #endif   /// 更新 row /= 2; col /=2; t++; src = dst; }   return dst; }   /// 小波逆变换 Mat IWDT( const Mat &_src, const string _wname, const int _level )const { int reValue = THID_ERR_NONE; Mat src = Mat_<float>(_src); Mat dst = Mat::zeros( src.rows, src.cols, src.type() ); int N = src.rows; int D = src.cols;   /// 高通低通滤波器 Mat lowFilter; Mat highFilter; wavelet( _wname, lowFilter, highFilter );   /// 小波变换 int t=1; int row = N/std::pow( 2., _level-1); int col = D/std::pow(2., _level-1);   while ( row<=N && col<=D ) { /// 小波列逆变换 for ( int j=0; j<col; j++ ) { /// 取出src数据的一行输入 Mat oneCol = Mat::zeros( row, 1, src.type() ); for ( int i=0; i<row; i++ ) { oneCol.at<float>(i,0) = src.at<float>(i,j); } oneCol = ( waveletReconstruct( oneCol.t(), lowFilter, highFilter ) ).t();   for ( int i=0; i<row; i++ ) { dst.at<float>(i,j) = oneCol.at<float>(i,0); } }   #if 0 //normalize( dst, dst, 0, 255, NORM_MINMAX ); IplImage dstImg2 = IplImage(dst); cvSaveImage( "dst.jpg", &dstImg2 ); #endif ///行小波逆变换 for( int i=0; i<row; i++ ) { /// 取出src中要处理的数据的一行 Mat oneRow = Mat::zeros( 1,col, src.type() ); for ( int j=0; j<col; j++ ) { oneRow.at<float>(0,j) = dst.at<float>(i,j); } oneRow = waveletReconstruct( oneRow, lowFilter, highFilter ); /// 将src这一行置为oneRow中的数据 for ( int j=0; j<col; j++ ) { dst.at<float>(i,j) = oneRow.at<float>(0,j); } }   #if 0 //normalize( dst, dst, 0, 255, NORM_MINMAX ); IplImage dstImg1 = IplImage(dst); cvSaveImage( "dst.jpg", &dstImg1 ); #endif   row *= 2; col *= 2; src = dst; }   return dst; }       /// 调用函数   /// 生成不同类型的小波，现在只有haar，sym2 void wavelet( const string _wname, Mat &_lowFilter, Mat &_highFilter )const { if ( _wname=="haar" || _wname=="db1" ) { int N = 2; _lowFilter = Mat::zeros( 1, N, CV_32F ); _highFilter = Mat::zeros( 1, N, CV_32F ); _lowFilter.at<float>(0, 0) = 1/sqrtf(N); _lowFilter.at<float>(0, 1) = 1/sqrtf(N);   _highFilter.at<float>(0, 0) = -1/sqrtf(N); _highFilter.at<float>(0, 1) = 1/sqrtf(N); } if ( _wname =="sym2" ) { int N = 4; float h[] = {-0.483, 0.836, -0.224, -0.129 }; float l[] = {-0.129, 0.224, 0.837, 0.483 };   _lowFilter = Mat::zeros( 1, N, CV_32F ); _highFilter = Mat::zeros( 1, N, CV_32F );   for ( int i=0; i<N; i++ ) { _lowFilter.at<float>(0, i) = l[i]; _highFilter.at<float>(0, i) = h[i]; }   } }   /// 小波分解 Mat waveletDecompose( const Mat &_src, const Mat &_lowFilter, const Mat &_highFilter )const { assert( _src.rows==1 && _lowFilter.rows==1 && _highFilter.rows==1 ); assert( _src.cols>=_lowFilter.cols && _src.cols>=_highFilter.cols ); Mat &src = Mat_<float>(_src);   int D = src.cols; Mat &lowFilter = Mat_<float>(_lowFilter); Mat &highFilter = Mat_<float>(_highFilter);     /// 频域滤波，或时域卷积；ifft( fft(x) * fft(filter)) = cov(x,filter) Mat dst1 = Mat::zeros( 1, D, src.type() ); Mat dst2 = Mat::zeros( 1, D, src.type() );   filter2D( src, dst1, -1, lowFilter ); filter2D( src, dst2, -1, highFilter );     /// 下采样 Mat downDst1 = Mat::zeros( 1, D/2, src.type() ); Mat downDst2 = Mat::zeros( 1, D/2, src.type() );   resize( dst1, downDst1, downDst1.size() ); resize( dst2, downDst2, downDst2.size() );     /// 数据拼接 for ( int i=0; i<D/2; i++ ) { src.at<float>(0, i) = downDst1.at<float>( 0, i ); src.at<float>(0, i+D/2) = downDst2.at<float>( 0, i ); }   return src; }   /// 小波重建 Mat waveletReconstruct( const Mat &_src, const Mat &_lowFilter, const Mat &_highFilter )const { assert( _src.rows==1 && _lowFilter.rows==1 && _highFilter.rows==1 ); assert( _src.cols>=_lowFilter.cols && _src.cols>=_highFilter.cols ); Mat &src = Mat_<float>(_src);   int D = src.cols;   Mat &lowFilter = Mat_<float>(_lowFilter); Mat &highFilter = Mat_<float>(_highFilter);   /// 插值; Mat Up1 = Mat::zeros( 1, D, src.type() ); Mat Up2 = Mat::zeros( 1, D, src.type() );   /// 插值为0 //for ( int i=0, cnt=1; i<D/2; i++,cnt+=2 ) //{ // Up1.at<float>( 0, cnt ) = src.at<float>( 0, i ); ///< 前一半 // Up2.at<float>( 0, cnt ) = src.at<float>( 0, i+D/2 ); ///< 后一半 //}   /// 线性插值 Mat roi1( src, Rect(0, 0, D/2, 1) ); Mat roi2( src, Rect(D/2, 0, D/2, 1) ); resize( roi1, Up1, Up1.size(), 0, 0, INTER_CUBIC ); resize( roi2, Up2, Up2.size(), 0, 0, INTER_CUBIC );   /// 前一半低通，后一半高通 Mat dst1 = Mat::zeros( 1, D, src.type() ); Mat dst2= Mat::zeros( 1, D, src.type() ); filter2D( Up1, dst1, -1, lowFilter ); filter2D( Up2, dst2, -1, highFilter );   /// 结果相加 dst1 = dst1 + dst2;   return dst1;   } =========================================================================== ***************************************************************************************************************** 其他代码实现 ***************************************************************************************************************** // 哈尔小波.cpp : 定义控制台应用程序的入口点。 // #include “stdafx.h” #include <iostream> #include “cv.h” #include “highgui.h” #include <ctype.h> // 二维离散小波变换（单通道浮点图像） void DWT(IplImage *pImage, int nLayer) {    // 执行条件    if (pImage)    {       if (pImage->nChannels == 1 &&          pImage->depth == IPL_DEPTH_32F &&          ((pImage->width >> nLayer) << nLayer) == pImage->width &&          ((pImage->height >> nLayer) << nLayer) == pImage->height)       {          int     i, x, y, n;          float   fValue   = 0;          float   fRadius  = sqrt(2.0f);          int     nWidth   = pImage->width;          int     nHeight  = pImage->height;          int     nHalfW   = nWidth / 2;          int     nHalfH   = nHeight / 2;          float **pData    = new float*[pImage->height];          float  *pRow     = new float[pImage->width];          float  *pColumn  = new float[pImage->height];          for (i = 0; i < pImage->height; i++)          {             pData[i] = (float*) (pImage->imageData + pImage->widthStep * i);          }          // 多层小波变换          for (n = 0; n < nLayer; n++, nWidth /= 2, nHeight /= 2, nHalfW /= 2, nHalfH /= 2)          {             // 水平变换             for (y = 0; y < nHeight; y++)             {                // 奇偶分离                memcpy(pRow, pData[y], sizeof(float) * nWidth);                for (i = 0; i < nHalfW; i++)                {                   x = i * 2;                   pData[y][i] = pRow[x];                   pData[y][nHalfW + i] = pRow[x + 1];                }                // 提升小波变换                for (i = 0; i < nHalfW - 1; i++)                {                   fValue = (pData[y][i] + pData[y][i + 1]) / 2;                   pData[y][nHalfW + i] -= fValue;                }                fValue = (pData[y][nHalfW - 1] + pData[y][nHalfW - 2]) / 2;                pData[y][nWidth - 1] -= fValue;                fValue = (pData[y][nHalfW] + pData[y][nHalfW + 1]) / 4;                pData[y][0] += fValue;                for (i = 1; i < nHalfW; i++)                {                   fValue = (pData[y][nHalfW + i] + pData[y][nHalfW + i - 1]) / 4;                   pData[y][i] += fValue;                }                // 频带系数                for (i = 0; i < nHalfW; i++)                {                   pData[y][i] *= fRadius;                   pData[y][nHalfW + i] /= fRadius;                }             }             // 垂直变换             for (x = 0; x < nWidth; x++)             {                // 奇偶分离                for (i = 0; i < nHalfH; i++)                {                   y = i * 2;                   pColumn[i] = pData[y][x];                   pColumn[nHalfH + i] = pData[y + 1][x];                }                for (i = 0; i < nHeight; i++)                {                   pData[i][x] = pColumn[i];                }                // 提升小波变换                for (i = 0; i < nHalfH - 1; i++)                {                   fValue = (pData[i][x] + pData[i + 1][x]) / 2;                   pData[nHalfH + i][x] -= fValue;                }                fValue = (pData[nHalfH - 1][x] + pData[nHalfH - 2][x]) / 2;                pData[nHeight - 1][x] -= fValue;                fValue = (pData[nHalfH][x] + pData[nHalfH + 1][x]) / 4;                pData[0][x] += fValue;                for (i = 1; i < nHalfH; i++)                {                   fValue = (pData[nHalfH + i][x] + pData[nHalfH + i - 1][x]) / 4;                   pData[i][x] += fValue;                }                // 频带系数                for (i = 0; i < nHalfH; i++)                {                   pData[i][x] *= fRadius;                   pData[nHalfH + i][x] /= fRadius;                }             }          }          delete[] pData;          delete[] pRow;          delete[] pColumn;       }    } } // 二维离散小波恢复（单通道浮点图像） void IDWT(IplImage *pImage, int nLayer) {    // 执行条件    if (pImage)    {       if (pImage->nChannels == 1 &&          pImage->depth == IPL_DEPTH_32F &&          ((pImage->width >> nLayer) << nLayer) == pImage->width &&          ((pImage->height >> nLayer) << nLayer) == pImage->height)       {          int     i, x, y, n;          float   fValue   = 0;          float   fRadius  = sqrt(2.0f);          int     nWidth   = pImage->width >> (nLayer - 1);          int     nHeight  = pImage->height >> (nLayer - 1);          int     nHalfW   = nWidth / 2;          int     nHalfH   = nHeight / 2;          float **pData    = new float*[pImage->height];          float  *pRow     = new float[pImage->width];          float  *pColumn  = new float[pImage->height];          for (i = 0; i < pImage->height; i++)          {             pData[i] = (float*) (pImage->imageData + pImage->widthStep * i);          }          // 多层小波恢复          for (n = 0; n < nLayer; n++, nWidth *= 2, nHeight *= 2, nHalfW *= 2, nHalfH *= 2)          {             // 垂直恢复             for (x = 0; x < nWidth; x++)             {                // 频带系数                for (i = 0; i < nHalfH; i++)                {                   pData[i][x] /= fRadius;                   pData[nHalfH + i][x] *= fRadius;                }                // 提升小波恢复                fValue = (pData[nHalfH][x] + pData[nHalfH + 1][x]) / 4;                pData[0][x] -= fValue;                for (i = 1; i < nHalfH; i++)                {                   fValue = (pData[nHalfH + i][x] + pData[nHalfH + i - 1][x]) / 4;                   pData[i][x] -= fValue;                }                for (i = 0; i < nHalfH - 1; i++)                {                   fValue = (pData[i][x] + pData[i + 1][x]) / 2;                   pData[nHalfH + i][x] += fValue;                }                fValue = (pData[nHalfH - 1][x] + pData[nHalfH - 2][x]) / 2;                pData[nHeight - 1][x] += fValue;                // 奇偶合并                for (i = 0; i < nHalfH; i++)                {                   y = i * 2;                   pColumn[y] = pData[i][x];                   pColumn[y + 1] = pData[nHalfH + i][x];                }                for (i = 0; i < nHeight; i++)                {                   pData[i][x] = pColumn[i];                }             }             // 水平恢复             for (y = 0; y < nHeight; y++)             {                // 频带系数                for (i = 0; i < nHalfW; i++)                {                   pData[y][i] /= fRadius;                   pData[y][nHalfW + i] *= fRadius;                }                // 提升小波恢复                fValue = (pData[y][nHalfW] + pData[y][nHalfW + 1]) / 4;                pData[y][0] -= fValue;                for (i = 1; i < nHalfW; i++)                {                   fValue = (pData[y][nHalfW + i] + pData[y][nHalfW + i - 1]) / 4;                   pData[y][i] -= fValue;                }                for (i = 0; i < nHalfW - 1; i++)                {                   fValue = (pData[y][i] + pData[y][i + 1]) / 2;                   pData[y][nHalfW + i] += fValue;                }                fValue = (pData[y][nHalfW - 1] + pData[y][nHalfW - 2]) / 2;                pData[y][nWidth - 1] += fValue;                // 奇偶合并                for (i = 0; i < nHalfW; i++)                {                   x = i * 2;                   pRow[x] = pData[y][i];                   pRow[x + 1] = pData[y][nHalfW + i];                }                memcpy(pData[y], pRow, sizeof(float) * nWidth);             }          }          delete[] pData;          delete[] pRow;          delete[] pColumn;       }    } } int _tmain(int argc, _TCHAR* argv[]) { // 小波变换层数 int nLayer = 1; // 输入彩色图像 IplImage *pSrc = cvLoadImage(“lena.bmp”, CV_LOAD_IMAGE_COLOR); // 计算小波图象大小 CvSize size = cvGetSize(pSrc); if ((pSrc->width >> nLayer) << nLayer != pSrc->width) {    size.width = ((pSrc->width >> nLayer) + 1) << nLayer; } if ((pSrc->height >> nLayer) << nLayer != pSrc->height) {    size.height = ((pSrc->height >> nLayer) + 1) << nLayer; } // 创建小波图象 IplImage *pWavelet = cvCreateImage(size, IPL_DEPTH_32F, pSrc->nChannels); if (pWavelet) {    // 小波图象赋值    cvSetImageROI(pWavelet, cvRect(0, 0, pSrc->width, pSrc->height));    cvConvertScale(pSrc, pWavelet, 1, -128);    cvResetImageROI(pWavelet);    // 彩色图像小波变换    IplImage *pImage = cvCreateImage(cvGetSize(pWavelet), IPL_DEPTH_32F, 1);    if (pImage)    {       for (int i = 1; i <= pWavelet->nChannels; i++)       {          cvSetImageCOI(pWavelet, i);          cvCopy(pWavelet, pImage, NULL);          // 二维离散小波变换          DWT(pImage, nLayer);          // 二维离散小波恢复          // IDWT(pImage, nLayer);          cvCopy(pImage, pWavelet, NULL);       }       cvSetImageCOI(pWavelet, 0);       cvReleaseImage(&pImage);    }    // 小波变换图象    cvSetImageROI(pWavelet, cvRect(0, 0, pSrc->width, pSrc->height));    cvConvertScale(pWavelet, pSrc, 1, 128);         cvNamedWindow(“pWavelet”,CV_WINDOW_AUTOSIZE);   cvShowImage(“pWavelet”,pWavelet);     cvNamedWindow(“pSrc”,CV_WINDOW_AUTOSIZE);   cvShowImage(“pSrc”,pSrc);   cvWaitKey(0);    cvResetImageROI(pWavelet); // 本行代码有点多余，但有利用养成良好的编程习惯    cvReleaseImage(&pWavelet); } // 显示图像pSrc // … cvReleaseImage(&pSrc); return 0; } ************************ ————————————-Codeproject Haar ————————————————————————————————- void Haar::transrows(char** dest, char** sour, unsigned int w, unsigned int h) const {         unsigned int w2 = w / 2;         __m64 m00FF;         m00FF.m64_u64 = 0x00FF00FF00FF00FF;         for (unsigned int y = 0; y < h; y++) {                 __m64 *mlo = (__m64 *) & dest[y][0];                 __m64 *mhi = (__m64 *) & dest[y][w2];                 __m64 *msour = (__m64 *) & sour[y][0];                 for (unsigned int k = 0; k < w2 / 8; k++) {   //k<w2/8   k=8*k                         __m64 even = _mm_packs_pu16(_mm_and_si64(*msour, m00FF), _mm_and_si64(*(msour + 1), m00FF));       //even coeffs                         __m64 odd = _mm_packs_pu16(_mm_srli_pi16(*msour, 8), _mm_srli_pi16(*(msour + 1), 8));              //odd coeffs                         addsub(even, odd, mlo++, mhi++);                         msour += 2;                 }                 if (w2 % 8) {                         for (unsigned int k = w2 - (w2 % 8); k < w2; k++) {                                 dest[y][k] = char(((int)sour[y][2*k] + (int)sour[y][2*k+1]) / 2);                                 dest[y][k+w2] = char(((int)sour[y][2*k] - (int)sour[y][2*k+1]) / 2);                         }                 }         }         _mm_empty(); } void Haar::transcols(char** dest, char** sour, unsigned int w, unsigned int h) const {         unsigned int h2 = h / 2;         for (unsigned int k = 0; k < h2; k++) {                 __m64 *mlo = (__m64 *) & dest[k][0];                 __m64 *mhi = (__m64 *) & dest[k+h2][0];                 __m64 *even = (__m64 *) & sour[2*k][0];                 __m64 *odd = (__m64 *) & sour[2*k+1][0];                 for (unsigned int x = 0; x < w / 8; x++) {                                                  addsub(*even, *odd, mlo, mhi);                                                 even++;                         odd++;                         mlo++;                         mhi++;                 }         }         _mm_empty();         //odd remainder         for (unsigned int x = w - (w % 8); x < w; x++) {                 for (unsigned int k = 0; k < h2; k++) {                                                 dest[k][x] = char(((int)sour[2*k][x] + (int)sour[2*k+1][x]) / 2);                                                 dest[k+h2][x] = char(((int)sour[2*k][x] - (int)sour[2*k+1][x]) / 2);                 }         } } **************************************************************************************************************** 哈尔小波变换示例 1.哈尔基函数   最简单的基函数是哈尔基函数(Haar basis function)。哈尔基函数在1909年提出，它是由一组分段常值函数组成的函数集。这个函数集定义在半开区间[0,1)上，每一个分段常值函数的数 值在一个小范围里是1，其他地方为0，现以图像为例并使用线性代数中的矢量空间来说明哈尔基函数。   如果一幅图像仅由2^0=1个像素组成，这幅图像在整个[0,1) 区间中就是一个常值函数。 用q_00(x)表示这个常值函数，用V0表示由这个常值函数生成的矢量空间， 即V0：q_00(x)=1(0<=x<1)或0(其它),它是构成矢量空间V0的基。   如果一幅图像由2^1=2个像素组成，这幅图像在[0,1) 区间中有两个等间隔的子区间：[0,1/2) 和[1/2,1) ，每一个区间中各有1个常值函数，分别用q_10(x) q_11(x)表示。用V1表示由2个子区间中的常值函数生成的矢量空间，即V1: q_10(x)  =1(0<=x<1/2)或0(其它); q_11(x)=1(1/2<=x<1)或0(其它).这2个常值函数就是构成矢量空间V1的基。   如果一幅图像由2^2= 4个像素组成，这幅图像在[0,1)区间中被分成4个等间隔的子区间：[0,1/4),[1/4,1/2),[1/2,3/4)和[3/4,1)，它们的 常值函数分别用q_20(x) q_21(x) q_22(x) q_23(x)表示，用V2表示由4个子区间中的常值函数生成的矢量空间，即V2:q_20(x)=1(0<=x<1/4)或0(其它); q_21(x)=1(1/4<=x<1/2)或0(其它); q_22(x) =1(1/2<=x<3/4)或0(其它); q_23(x)=1(3/4<=x<1)或0(其它).这4个常值函数就是构成矢量空间V2的基。   我们可以按照这种方法继续定义基函数和由它生成的矢量空间。总而言之，为了表示矢量空间中的矢量，每一个矢量空间V_i都需要定义一个基(basis)。 为生成矢量空间V_i而定义的基函数也叫做尺度函数,这种函数通常用符号q_ij(x)表示。哈尔基函数定义为q(x)=1(0<=x<1) 或0(其它),哈尔尺度函数q_ij(x)定义为q_ij(x)=q(2^i·x-j),j=0,1,…, 2^i-1其中,i为尺度因子，改变i使函数图形缩小或者放大,j为平移参数，改变j使函数沿x轴方向平移。   空间矢量V_i定义为V_i=span{ q_ij(x)},j=0,1,…, 2^i-1.由于定义了基和矢量空间，我们就可以把由2^i个像素组成的一维图像看成为矢量空间V_i中的矢量。由于这些矢量都是在单位区间[0,1)上 定义的函数，所以在V_i矢量空间中的每一个矢量也被包含在V_i+1矢量空间中，即V0∈V1∈…∈V_i∈V_i+1. 2.哈尔小波函数   小波函数通常用Ψ_ij(x)表示。与哈尔基函数相对应的小波称为基本哈尔小波函数,并由下式定义: Ψ(x)=1(0<=x<1/2)或 -1(1/2<=x<1)或0(其它);  哈尔小波尺度函数Ψ_ij(x)定义为 Ψ_ij(x)= Ψ(2^i·x-j),j=0,1,…, 2^i-1.用小波函数构成的矢量空间Wi表示，Wi=span{Ψ_ij(x)}, j=0,1,…, 2^i-1   根据哈尔小波函数的定义，可以写出生成W0，W1和W2等矢量空间的小波函数。 生成矢量空间W0的哈尔小波：Ψ_00(x)=1(0<=x<1/2)或-1(1/2<=x<1)或0(其它) 生成矢量空间W1的哈尔小波：Ψ_10(x)=1(0<=x<1/4)或-1(1/4<=x<1/2)或0(其它); Ψ_11(x)=1(1/2<=x<3/4)或-1(3/4<=x<1/2)或0(其它). 生成矢量空间W2的哈尔小波：Ψ_00(x)=1(0<=x<1/2)或-1(1/2<=x<1)或0(其它) Ψ_20(x)=1(0<=x<1/8)或-1(1/8<=x<2/8)或0(其它); Ψ_21(x)=1(2/8<=x<3/8)或-1(3/8<=x<4/8)或0(其它); Ψ_22(x)=1(4/8<=x<5/8)或-1(5/8<=x<6/8)或0(其它); Ψ_00(x)=1(6/8<=x<7/8)或-1(7/8<=x<1)或0(其它). 哈尔基和哈尔小波分别使用下面两个式子进行规范化 q_ij(x)=2^(i/2)·q(2^i·x-j), Ψ_ij(x)= 2^(i/2)·Ψ(2^i·x-j), 3.哈尔基的结构   使用哈尔基函数 q_ij(x)和哈尔小波函数Ψ_ij(x)生成的矢量空间Vi和Wi具有下面的性质，V_i+1=V_i⊕W_i。这就是说，在矢量空间V_i+1中, 矢量空间W_i中的小波可用来表示一个函数在矢量空间V_i 中不能表示的部分。因此，小波可定义为生成矢量空间W_i的一组线性无关的函数Ψ_ij(x)的集合。 4.哈尔小波变换   小波变换的基本思想是用一组小波函数或者基函数表示一个函数或者信号，例如图像信号。为了理解什么是小波变换，下面用一个具体的例子来说明小波变换的过程。假设有一幅分辨率只有4个像素的一维图像，对应的像素值分别为[9 7 3 5]   计算它的哈尔小波变换系数: (1).求均值(averaging)。计算相邻像素对的平均值，得到一幅分辨率比较低的新图像，它的像素数目变成了2个，即新的图像的分辨率是原来的1/2，相应的像素值为：[8 4] (2). 求差值(diqqerencing)。很明显，用2个像素表示这幅图像时，图像的信息已经部分丢失。为了能够从由2个像素组成的图像重构出由4个像素组成 的原始图像，就需要存储一些图像的细节系数,以便在重构时找回丢失的信息。方法是把像素对的第一个像素值减去这个像素对的平均值，或者使用这个像素对的差 值除以2。在这个例子中，第一个细节系数是(9-8)=1，因为计算得到的平均值是8,它比9小1而比7大1，存储这个细节系数就可以恢复原始图像的前两 个像素值。使用同样的方法，第二个细节系数是(3-4)=-1，存储这个细节系数就可以恢复后2个像素值。因此，原始图像就可以用下面的两个平均值和两个 细节系数表示[8 4 1 -1] (3).重复第1，2步，把由第一步分解得到的图像进一步分解成分辨率更低的图像[6 2 1 -1].   由此可见，通过上述分解就把由4像素组成的一幅图像用一个平均像素值和三个细节系数表示，这个过程就叫做哈尔小波变换,也称哈尔小波分解(Haar wavelet decomposition)。这个概念可以推广到使用其他小波基的变换。在这个例子中我们可以看到，①变换过程中没有丢失信息，因为能够从所记录的数据 中重构出原始图像。②对这个给定的变换，我们可以从所记录的数据中重构出各种分辨率的图像。例如，在分辨率为1的图像基础上重构出分辨率为2的图像，在分 辨率为2的图像基础上重构出分辨率为4的图像。③通过变换之后产生的细节系数的幅度值比较小，这就为图像压缩提供了一种途径，例如去掉一些微不足道的细节 系数而不影响对重构图像的理解。   求均值和差值的过程实际上就是一维小波变换的过程，现在用数学方法重新描述小波变换的过程。 (1) 用V2中的哈尔基表示 图像I(x)=[9 7 3 5]有2^j =2^2=4像素，因此可以用生成矢量空间V2中的基函数的线性组合表示， I(x) =9 q_20(x)+ 7q_21(x)+3q_22(x)+5q_23(x) (2)用V1和W1中的函数表示 根据V2=V1⊕W1,I(x)可表示成 I(x)=8q_10(x)+4q_11(x)+Ψ_10(x)-Ψ_11(x) (3)用V0,W0和W1中的函数表示 V2=V0⊕W0⊕W1,I(x)可表示成 I(x) =6 q_00(x)+2Ψ_00(x)+Ψ_10(x)-Ψ_11(x) 5.规范化算法   规范化的小波变换与非规范化的小波变换相比，唯一的差别是在规范化的小波变换中需要乘一个规范化的系数。下面用一个例子说明。对函数f(x)= [2,5,8,9,7,4, -1, 1]作哈尔小波变换。哈尔小波变换实际上是使用求均值和差值的方法进行分解。我们把f(x)看成是矢量空间V3中的一个向量，尺度因子i= 3，因此最多可分解为3层.分解过程如下。 步骤1：c=(2 5,8 9,7,4,-1,1)/√2=(7,17,11,0, -3, -1,3, -2)/√2 解释： 平均值（7,17,11,0）/2  差值（-3，-1,3，-2）/2（与平均值的差值相对于直接差值的1/2  即 A- (A+B)/2 = (A-B)/2  ,我们将除以2提出来） 所以步骤一的结果是（7,17,11,0）/2  并上（-3，-1,3，-2）/2 = (7,17,11,0, -3, -1,3, -2)/2 步骤2：c=[(7+17)/√2,(11+0)/√2,(7-17)/√2,(11-0)/√2,-3,-1,3,-2]/√2   =(24/√2,11/√2,-10/√2,11/√2,-3,-1,3,-2)/√2 步骤3:c=[(24+11)/2,(24-11)/2 ,-10/√2,11/√2,-3,-1,3,-2]/√2    =(35/2,13/2,-10/√2,11/√2,-3,-1,3,-2)/√2 这个结果实际上是后面哈尔转换N=8转换矩阵相乘的结果。 N = 8 ，