离散傅立叶变换(DFT)

 void cvShiftDFT(CvArr * src_arr, CvArr * dst_arr )
{
    CvMat * tmp;
    CvMat q1stub, q2stub;
    CvMat q3stub, q4stub;
    CvMat d1stub, d2stub;
    CvMat d3stub, d4stub;
    CvMat * q1, * q2, * q3, * q4;
    CvMat * d1, * d2, * d3, * d4;

    CvSize size = cvGetSize(src_arr);
    CvSize dst_size = cvGetSize(dst_arr);
    int cx, cy;

    if(dst_size.width != size.width ||
       dst_size.height != size.height){
        cvError( CV_StsUnmatchedSizes, "cvShiftDFT", "Source and Destination arrays must have equal sizes", __FILE__, __LINE__ );  
    }

    if(src_arr==dst_arr){
        tmp = cvCreateMat(size.height/2, size.width/2, cvGetElemType(src_arr));
    }
   
    cx = size.width/2;
    cy = size.height/2; // image center

    q1 = cvGetSubRect( src_arr, &q1stub, cvRect(0,0,cx, cy) );
    q2 = cvGetSubRect( src_arr, &q2stub, cvRect(cx,0,cx,cy) );
    q3 = cvGetSubRect( src_arr, &q3stub, cvRect(cx,cy,cx,cy) );
    q4 = cvGetSubRect( src_arr, &q4stub, cvRect(0,cy,cx,cy) );
    d1 = cvGetSubRect( src_arr, &d1stub, cvRect(0,0,cx,cy) );
    d2 = cvGetSubRect( src_arr, &d2stub, cvRect(cx,0,cx,cy) );
    d3 = cvGetSubRect( src_arr, &d3stub, cvRect(cx,cy,cx,cy) );
    d4 = cvGetSubRect( src_arr, &d4stub, cvRect(0,cy,cx,cy) );

    if(src_arr!=dst_arr){
        if( !CV_ARE_TYPES_EQ( q1, d1 )){
            cvError( CV_StsUnmatchedFormats, "cvShiftDFT", "Source and Destination arrays must have the same format", __FILE__, __LINE__ );
        }
        cvCopy(q3, d1, 0);
        cvCopy(q4, d2, 0);
        cvCopy(q1, d3, 0);
        cvCopy(q2, d4, 0);
    }
    else{
        cvCopy(q3, tmp, 0);
        cvCopy(q1, q3, 0);
        cvCopy(tmp, q1, 0);
        cvCopy(q4, tmp, 0);
        cvCopy(q2, q4, 0);
        cvCopy(tmp, q2, 0);
    }
}

void change()
{
    IplImage * im;

    IplImage * realInput;
    IplImage * imaginaryInput;
    IplImage * complexInput;
    int dft_M, dft_N;
    CvMat* dft_A, tmp;
    IplImage * image_Re;
    IplImage * image_Im;
    double m, M;

    im = cvLoadImage("first2.bmp", CV_LOAD_IMAGE_GRAYSCALE );
    if( !im )
        return;

    realInput = cvCreateImage( cvGetSize(im), IPL_DEPTH_64F, 1);
    imaginaryInput = cvCreateImage( cvGetSize(im), IPL_DEPTH_64F, 1);
    complexInput = cvCreateImage( cvGetSize(im), IPL_DEPTH_64F, 2);

    cvScale(im, realInput, 1.0, 0.0);
    cvZero(imaginaryInput);
    cvMerge(realInput, imaginaryInput, NULL, NULL, complexInput);

    dft_M = cvGetOptimalDFTSize( im->height - 1 );
    dft_N = cvGetOptimalDFTSize( im->width - 1 );

    dft_A = cvCreateMat( dft_M, dft_N, CV_64FC2 );
    image_Re = cvCreateImage( cvSize(dft_N, dft_M), IPL_DEPTH_64F, 1);
    image_Im = cvCreateImage( cvSize(dft_N, dft_M), IPL_DEPTH_64F, 1);

    // copy A to dft_A and pad dft_A with zeros
    cvGetSubRect( dft_A, &tmp, cvRect(0,0, im->width, im->height));
    cvCopy( complexInput, &tmp, NULL );
    if( dft_A->cols > im->width )
    {
        cvGetSubRect( dft_A, &tmp, cvRect(im->width,0, dft_A->cols - im->width, im->height));
        cvZero( &tmp );
    }

    // no need to pad bottom part of dft_A with zeros because of
    // use nonzero_rows parameter in cvDFT() call below

    cvDFT( dft_A, dft_A, CV_DXT_FORWARD, complexInput->height );

    cvNamedWindow("win", 0);
    cvNamedWindow("magnitude", 0);
    cvShowImage("win", im);

    // Split Fourier in real and imaginary parts
    cvSplit( dft_A, image_Re, image_Im, 0, 0 );

    // Compute the magnitude of the spectrum Mag = sqrt(Re^2 + Im^2)
    cvPow( image_Re, image_Re, 2.0);
    cvPow( image_Im, image_Im, 2.0);
    cvAdd( image_Re, image_Im, image_Re, NULL);
    cvPow( image_Re, image_Re, 0.5 );

    // Compute log(1 + Mag)
    cvAddS( image_Re, cvScalarAll(1.0), image_Re, NULL ); // 1 + Mag
    cvLog( image_Re, image_Re ); // log(1 + Mag)


    // Rearrange the quadrants of Fourier image so that the origin is at
    // the image center
    cvShiftDFT( image_Re, image_Re );

    cvMinMaxLoc(image_Re, &m, &M, NULL, NULL, NULL);
    cvScale(image_Re, image_Re, 1.0/(M-m), 1.0*(-m)/(M-m));
    cvShowImage("magnitude", image_Re);

    cvWaitKey(-1);
}

### 离散傅里叶变换DFT原理 离散傅里叶变换(DFT)是一种用于分析有限长度序列的方法,它能够将时域中的信号转换成频域表示形式。对于一个长度为 \( N \)离散时间信号 \( x[n] \),其对应的 DFT 定义如下: \[ X[k]=\sum_{n=0}^{N-1}{x[n]\cdot e^{-j2\pi kn/N}} , k = 0, ..., N-1 \] 这里 \( j=\sqrt{-1} \), 表达式中指数项代表旋转因子。 通过上述公式可以计算得到输入序列的各个频率成分及其幅值大小[^1]。 当执行逆离散傅里叶变换(IDFT)时,则是从频域恢复原始的时间序列: \[ x[n]=(1/N)\sum_{k=0}^{N-1}{X[k]\cdot e^{j2\pi nk/N}}, n = 0,...,N−1 \] 这表明傅里叶变换本质上实现了数据从时域到频域之间的相互转化过程[^2]。 另外值得注意的是,在实际操作过程中为了更直观地观察频谱分布情况,通常会做一次频谱平移使得中心位于零频率处即[-π, π]区间内。 ### 应用实例 在工程实践中,DFT有着广泛的应用场景之一便是滤波器设计以及音频处理等领域。比如去除噪声干扰、提取特定声音特征等任务都可以借助该工具完成。具体来说,通过对含有杂音的声音文件实施快速傅立叶变化(FFT算法实现高效版DFT),分离出有用的信息部分后再经由相应的反向运算重建纯净版本。 此外,图像压缩技术JPEG标准也采用了类似的思路来减少冗余信息量从而达到节省存储空间的目的。在此基础上发展起来的小波变换更是进一步拓展了此类方法论的应用范围[^3]。 ```python import numpy as np from matplotlib import pyplot as plt # 创建测试信号 fs = 1000 # Sampling frequency t = np.linspace(0, 1, fs, endpoint=False) frequencies = [50, 120] amplitudes = [0.7, 1] signal = sum([amp * np.sin(2*np.pi*freq*t) for freq, amp in zip(frequencies, amplitudes)]) # 计算并绘制DFT结果 fft_result = np.fft.fft(signal)/len(t) frequency_domain = np.fft.fftfreq(len(t), d=1/fs) plt.figure(figsize=(8,6)) plt.plot(frequency_domain[:int(fs/2)], abs(fft_result)[:int(fs/2)]) plt.title('Magnitude Spectrum') plt.xlabel('Frequency(Hz)') plt.ylabel('|Amplitude|') plt.grid(True) plt.show() ``` 此段代码展示了如何利用 Python 中 `numpy` 和 `matplotlib` 库来进行简单的正弦波合成实验,并对其应用 FFT 来获取频谱图展示效果。
评论 1
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值