Windows LiveWindows Live™
 
 
Kailiang's profileoatmeal notesPhotosBlogListsMore Tools Help

Blog
    June 19

    fourier transform and MATLAB

    我不在摊子上的时候，他们竟将我的喷雾哩水和染发剂卖给了一个大爷。大爷可能头发会变白，需要染发，还问过染发剂什么颜色，他们说是紫黑色，实际上是紫红色，年轻人染比较合适。他们把六块钱给我的时候，觉得有些过意不去。不过人已经找不到了，只能作罢。中午很宏观，两个百吉馍，还有三鲜米线，黄瓜条，酸梅汁。罪过罪过，以后要少食多餐，才能陶冶情操，身体健康。晚上巴西进了两个球，又精神了，试了试傅立叶变换，发现这个东西挺好用，把信号变到频域，完了还能变回去。还发现一个很有用的例子，贴在下面。
    ----------------------------------------------------------------------------------------------------

    A common use of Fourier transforms is to find the frequency components of a signal buried in a noisy time domain signal. Consider data sampled at 1000 Hz. Form a signal containing 50 Hz and 120 Hz and corrupt it with some zero-mean random noise:

    • t = 0:0.001:0.6;
x = sin(2*pi*50*t)+sin(2*pi*120*t);
y = x + 2*randn(size(t));
plot(1000*t(1:50),y(1:50))
title('Signal Corrupted with Zero-Mean Random Noise')
xlabel('time (milliseconds)')

      

       

    It is difficult to identify the frequency components by looking at the original signal. Converting to the frequency domain, the discrete Fourier transform of the noisy signal y is found by taking the 512-point fast Fourier transform (FFT):

    • Y = fft(y,512);

    The power spectrum, a measurement of the power at various frequencies, is

    • Pyy = Y.* conj(Y) / 512;

    Graph the first 257 points (the other 255 points are redundant) on a meaningful frequency axis: 

    • f = 1000*(0:256)/512;
plot(f,Pyy(1:257))
title('Frequency content of y')
xlabel('frequency (Hz)')

      • ------------------------------------------------------------------------

      本来很规则的信号，加上一点噪声后，波形就完全看不出什么样子了， 不过通过傅立叶变换，就可以把原来信号的频率检出来。
      在MATLAB中，函数fft计算一个信号的离散傅立叶变换。在数据的长度是2的幂次或质因数的乘积的情况下，就用快速傅立叶变换（FFT)来计算离散傅立叶变换。当数据长度是2的幂次时，计算速度显著增加，因此，只要可能，选择数据长度为2的幂次或者用零来填补数据，使得数据长度等于2的幂次就非常重要。

      ------------------------------------------------------------------------

      Basic Operations                                                  
      ----------------------------------------------------------------
      cumprod    Cumulative product of elements.                 
      cumsum     Cumulative sum of elements.                     
      max        Largest component.                              
      mean       Average or mean value.                          
      median     Median value.                                   
      min        Smallest component.                             
      prod       Product of elements.                            
      sort       Sort in ascending order.                        
      std        Standard deviation.                             
      sum        Sum of elements.                                
      trapz      Numerical integration using trapezoidal method. 
      ----------------------------------------------------------------

      -------------------------------------------------------------------
      Finite Differences                                                  
      -------------------------------------------------------------------
      del2       Five-point discrete Laplacian.                  
      diff       Difference function and approximate derivative. 
      gradient   Approximate gradient.                           
      -------------------------------------------------------------------

      -------------------------------------
      Correlation                            
      -------------------------------------
      corrcoef   Correlation coefficients. 
      cov        Covariance matrix.        
      -------------------------------------

      ---------------------------------------------------------------------
      Filtering and Convolution                                             
      ---------------------------------------------------------------------
      conv       Convolution and polynomial multiplication. 
      conv2      Two-dimensional convolution.               
      deconv     Deconvolution and polynomial division.     
      filter     One-dimensional digital filter.            
      filter2    Two-dimensional digital filter.            
      ---------------------------------------------------------------------

      ----------------------------------------------------------------------
      Fourier Transforms                                                      
      ----------------------------------------------------------------------
      abs        Magnitude.                                          
      angle      Phase angle.                                        
      cplxpair   Sort numbers into complex conjugate pairs.          
      fft        Discrete Fourier transform.                         
      fft2       Two-dimensional discrete Fourier transform.         
      fftshift   Move zeroth lag to center of spectrum.              
      ifft       Inverse discrete Fourier transform.                 
      ifft2      Two-dimensional inverse discrete Fourier transform. 
      nextpow2   Next higher power of 2.                             
      unwrap     Remove phase angle jumps across 360\xa1  boundaries.    
      ----------------------------------------------------------------------

      ---------------------------------------
      Vector Functions                        
      ---------------------------------------
      cross      Vector cross product. 
      dot        Vector dot product.   
      ---------------------------------------

       

      精通MATLAB的傅立叶级数函数
       
      fsderiv(Kn,Wo)
       傅立叶级数的微分
       
      fseval(Kn,t,Wo)
       计算傅立叶级数
       
      fsfind(‘ fname ‘,T,N)
       寻找时间函数的傅立叶级数的系数
       
      [An,Bn,Ao]=fsform(Kn)

      Kn=fsform(An,Bn,Ao)
       傅立叶级数不同形式之间的转换
       
      fsharm(Kn,i)
       提取特殊的傅立叶级数的谐波
       
      fsmsv(Kn)
       计算信号的均方根值
       
      fsresize(Kn,N)
       重新调整傅立叶级数的系数向量的大小
       
      fsresp(Num,Den,Un,Wo)
       线性系统对输入傅立叶级数Un的富里哀级数响应
       
      fsround(Kn)
       设置次要的傅立叶级数的系数为0
       
      fswindow(N, ‘ type ‘)

      fswindow(Kn, ‘ type ‘)
       产生一个窗口的函数，使吉布（Gibb）现象极小

       

    3:09 AM | Blog it | 学习

    Comments (5)

    Please wait...
    Sorry, the comment you entered is too long. Please shorten it.
    You didn't enter anything. Please try again.
    Sorry, we can't add your comment right now. Please try again later.
    To add a comment, you need permission from your parent. Ask for permission
    Your parent has turned off comments.
    Sorry, we can't delete your comment right now. Please try again later.
    You've exceeded the maximum number of comments that can be left in one day. Please try again in 24 hours.
    Your account has had the ability to leave comments disabled because our systems indicate that you may be spamming other users. If you believe that your account has been disabled in error please contact Windows Live support.
    Complete the security check below to finish leaving your comment.
    The characters you type in the security check must match the characters in the picture or audio.

    To add a comment, sign in with your Windows Live ID (if you use Hotmail, Messenger, or Xbox LIVE, you have a Windows Live ID). Sign in


    Don't have a Windows Live ID? Sign up
    Kai Zhangwrote:
    玩光学必须搞傅立叶变换，低级是信号的主体，而高级是信号的细节。变换了，就什么都好分析了。
    July 7
    Leanne Liuwrote:
    Fourier Transform……前一阵子为它费了不少血~~~~
    原来你是搞这种东西的……
    July 2
    shanwrote:
    哎呀呀，羡慕，佩服
    June 26
    Erin Jiaowrote:
    以前学信号与系统的时候写的一个
    从调占空比和谐波次数可以看到吉布斯现象的说
     
    另外，小波变换与傅里叶变换zz
     

    如果有人问我，如果傅里叶变换没有学好（深入理解概念），是否能学好小波。答案是否定的。如果有人还问我，如果第一代小波变换没学好，能否学好第二代小波变换。答案依然是否定的。但若你问我，没学好傅里叶变换，能否操作（编程）小波变换，或是没学好第一代小波，能否操作二代小波变换，答案是肯定的。

    一、一、基的概念

    我们要明确的是基的概念。两者都是基，信号都可以分成无穷多个他们的和（叠加）。而展开系数就是基与信号之间的内积，更通俗的说是投影。展开系数大的，说明信号和基，是足够相似的。这也就是相似性检测的思想。但我们必须明确的是，傅里叶是0-2pi标准正交基，而小波是-inf到inf之间的基。因此，小波在实轴上是紧的。而傅里叶的基（正弦或余弦），与此相反。而小波能不能成为Reisz基，或标准稳定的正交基，还有其它的限制条件。此外，两者相似的还有就是PARSEVAL定理。（时频能量守恒）。

          

    二、二、离散化的处理

    傅里叶变换，是一种数学的精妙描述。但计算机实现，却是一步步把时域和频域离散化而来的。第一步，时域离散化，我们得到离散时间傅里叶变换（DTFT），频谱被周期化；第二步，再将频域离散化，我们得到离散周期傅里叶级数（DFS），时域进一步被周期化。第三步，考虑到周期离散化的时域和频域，我们只取一个周期研究，也就是众所周知的离散傅里叶变换（DFT）。这里说一句，DFT是没有物理意义的，它只是我们研究的需要。借此，计算机的处理才成为可能。

    下面我们谈谈小波。所有满足容许性条件（从-INF到+INF积分为零）的函数，都可以成为小波。小波作为尺度膨胀和空间移位的一组函数也就诞生了。但连续取值的尺度因子和平移因子，在时域计算量和频域的混叠来说，都是极为不便的。用更为专业的俗语，叫再生核。也就是，对于任何一个尺度a和平移因子b的小波，和原信号内积，所得到的小波系数，都可以表示成，在a，b附近生成的小波，投影后小波系数的线性组合。这就叫冗余性。这时的连续小波是与正交基毫无关系的东西，它顶多也只能作为一种积分变换或基。但它的显微镜特点和相似性检测能力，已经显现出来了。为了进一步更好的将连续小波变换离散化，以下步骤是一种有效方法。第一步，尺度离散化。一般只将a二进离散化，此时b是任意的。这样小波被称为二进小波。第二步，离散b。怎么离散化呢？b取多少才合适呢？于是，叫小波采样定理的东西，就这样诞生了。也就是小波平移的最小距离（采样间隔），应该大于二倍小波基的最高频率（好像类似，记不清了）。所以b取尺度的整数倍就行了。也就是越胖的小波，对应频谱越窄，平移量应该越大，采样间隔越大。当然，第一二两步的频域理解，即在满足频域窗口中心是3倍的频域窗口半径的前提下，频域就在统计上是完美二分的。（但很多小波满足不了这个条件，而且频域窗口能量不集中，所以只是近似二分的）。这时的小波变换，称为离散二进小波变换。第三步，引入稳定性条件。也就是经过变换后信号能量和原信号能量有什么不等式关系。满足稳定性条件后，也就是一个小波框架产生成了可能。他是数值稳定性的保证。一个稍弱的稳定条件，就是0<A<=B<+INF，并且小波函数线性无关，此时小波基称为Reisz基。并且，如果变换后能量守恒，(A=B=1)，并且线性无关，这就是标准离散正交小波基。这种分解也就是大家熟知的直和分解。若A和B不相等，且相差很大，我们就说小波不是紧框架的，所以双正交，对偶小波也就自然而然引进来了。若A和B不相等，但又相差不大，这时稳定重构也是可能的，这时成为几乎紧框架的。（好像说这样小波有橹棒性特点，也就是粗略分解，但却精确重构。）经过3步，我们最终地得到了一个二进离散化稳定的小波变换，这正是我们要的结果。

    三、三、快速算法

    如果说现代数字信号处理革命的算法，甚至是很多快速算法的老始祖，或者是满矩阵向量乘法一个几乎不可抗拒的最小计算量NlogN，那就是令我不得不佩服的快速傅里叶变换（FFT）。这里我不想解释过多的基2算法，和所谓的三重循环，还有那经典的蝶形单元，或是分裂基之类，我想说的就是一种时频对应关系。也就是算法的来源。我们首先明确，时域的卷积对应频域的相乘，因此我们为了实现卷积，可以先做傅里叶变换，接着在频域相乘，最后再做反傅里叶变换。这里要注意，实际我们在玩DSP。因此，大家要记住，圆周卷积和离散傅里叶变换，是一家子。快速傅里叶是离散傅里叶的快速算法。因此，我们实现离散线性卷积，先要补零。然后使得它和圆周卷积相等。然后就是快速傅里叶变换，频域相乘，最后反快速傅里叶变换。当然，如果我们就需要的是圆周卷积，那我们也就不需要多此一举的补零。这里，我们可以把圆周卷积，写成矩阵形式。这点很重要。Y=AX。这里的A是循环矩阵。但不幸的是A仍然是满阵。

    小波的快速算法。MALLET算法，是一个令人振奋的东西。它实质给了多分辨率分析（多尺度分析）一个变得一发而不可收的理由。它实质上，讲了这样一个意思。也就是。我在一个较高的尺度（细节）上作离散二进稳定的小波变换，得到了一个结果（小波系数），我若是想得到比它尺度低的小波系数（概貌），我不用再计算内积，只是把较高尺度的小波系数和低通或高通滤波器卷积再抽取即可。但是，所有这些证明的推导是在整个实轴上进行的。即把信号看成无限长的。但这仍不是我们想要的。还有，我们还必须在较高尺度上作一次内积，才可以使用此算法。因此，我们开始简化，并扩展此理论。第一，我们把信号的采样，作为一个较高层的小波系数近似初始值。（这是可以的，因为小波很瘦时，和取样函数无异）。第二，我们把原来的卷积，换为圆周卷积。这和DSP何尝不一样呢？他的物理意义，就是把信号作周期延拖（边界处理的一种），使之在整个实轴上扩展。这种算法令我为之一贯坚持的是，它是完全正交的，也就是说是正交变换。正变换Y=AX；反变换X=A’Y；一般对于标准正交基，A’是A的共轭转置，对于双正交A’是A的对偶矩阵。但不管如何，我们可以大胆的写，AA’=A’A=I。这里I是单位矩阵。

    那怎样操作才是最快的呢？我们来分析A的特点，首先A是正交阵，其次A是有循环矩阵特点，但此时A上半部分是由低通滤波器构成的循环子矩阵，下半部分是由高通滤波器构成的子矩阵，但却是以因子2为循环的。为什么，因为你做了2抽取。所以我们可以，实现小波变换用快速傅里叶变换。这时如果A是满阵的，则复杂度由O(N.^2）下降到O(NlogN)。（这个程序我已经传在了研学上，在原创区）。但还有一点，我们忘了A是稀疏的，因为信号是很长的，而滤波器确实很短的，也就是这个矩阵是个近似对角阵。所以，快速傅里叶是不快的，除非你傻到含有零的元素，也作了乘法。因此，小波变换是O(N)复杂度的。这是它的优势。但要实现，却不是那么容易，第一个方法，稀疏矩阵存储和稀疏矩阵乘法。第二个方法，因子化。因子化，是一个杰出的贡献。它在原有的O(N)的复杂度基础上，对于长滤波器，又把复杂度降低一半。但量级仍然是O(N)。

    四、四、时频分析

    对于平稳信号，傅里叶再好不过了。它反映的是信号总体的整个时间段的特点。在频率上，是点频的。而对于非平稳信号，它就无能为力了。而小波恰好对此派上用场。小波是反映信号，某个时间段的特点的。在频域上，是某个频率段的表现。但小波，作为频谱分析确实存在很多问题。但我们确实可以做出很多的小波满足这个特点。大家可以看冉启文的《小波变换与分数傅里叶变换》书，这里我不再赘述。还有，我们老是说小波是近似频域二分的，这在DSP上是怎样的，最近我也在思考。

    五、压缩、消噪、特征提取

           傅里叶变换的压缩，已经广泛应用了。它的简化版本就是DCT变换。而小波包的提出，也就使DCT有些相形见拙。首先，它提出代价函数，一般就是熵准则。其次，一个自适应树分解。再次，基于矩阵范数或较少位编码的稀疏化策略。这些使小波包的压缩近乎完美。小波包是从频域上实现的。从时域上，我们也可采用类似的分裂和并算法，来实现信号最优的表达，这种可变窗小波成为MALVAR小波。记住，压缩是小波最大的优势。

           消噪，一般的傅里叶算法，一般可以是IIR滤波和FIR滤波。两者各有优缺点。而小波的消噪，一般也是由多层分解和阈值策略组成。我们需要的是信号的特点，噪声的特点，然后确定用不用小波，或用什么小波。这点上，小波的优势并不是很明显。

           特征提取。这是小波的显微镜特点很好地运用。利用模极大值和LIPSCHITZ指数，我们可以对信号的突变点做分析。但这里面的问题也是很多。首先，在不同尺度上，噪声和信号的模极大值变化不同。再次，一般我们用求内积方法，求模极大值，而不用MALLET算法，或者改用叫多孔算法的东西来做。这点，我没任何体会，希望大家多讨论吧。

          这里，我不能谈应用很多的细节。但我们必须明确：1。你要对小波概念有着明确的理解。对诸如多分辨率，时频窗口与分析，框架，消失矩，模极大值，LIPSCHITZ指数等有着清醒地认识。2。你必须考虑小波在此问题上的可行性，这点尤为重要，小波不是万能的。

    June 19
    Erin Jiaowrote:
    function f=jjqff
    clc
    t=-1:0.001:3;
    a=input('请输入占空比:a=');
    n=input('请输入谐波次数:n=');
    f=(sign(t)-sign(t-a)+sign(t-1)-sign(t-1-a)+sign(t-2)-sign(t-2-a))/2;
    a0=a;
    for i=1:n
      
        an(1,i)=(sin(2*pi*i*a))/(pi*i);
        bn(1,i)=(1-cos(2*pi*i*a))/(pi*i);
    end
    fo=a0;
    for i=1:n
        F(1,i)=0.5*(an(1,i)^2+bn(1,i)^2)^(1/2);
        fi(1,i)=-atan(bn(1,i)/an(1,i));
        fo=an(1,i)*cos(i*2*pi*t)+bn(1,i)*sin(i*2*pi*t)+fo;
    end
    N=1:n;
    N=2*pi*N;
    subplot(2,2,1),plot(t,f),axis([-0.4,3,-0.4,1.4]),title('原波')
    subplot(2,2,3),plot(t,fo),axis([-0.4,3,-0.4,1.4]),title('谐波叠加')
    subplot(2,2,2),stem(N,F),title('单边频谱图')
    subplot(2,2,4),stem(N,fi),title('单边相位图')
    June 19

    Trackbacks

    The trackback URL for this entry is:
    http://oatmeal23.spaces.live.com/blog/cns!351CC67D458396D8!196.trak
    Weblogs that reference this entry
    • None