Python机器学习（五十六）SciPyfftpack（傅里叶变换）

　　SciPy提供了fftpack模块，包含了傅里叶变换的算法实现。
　　傅里叶变换把信号从时域变换到频域，以便对信号进行处理。傅里叶变换在信号与噪声处理、图像处理、音频信号处理等领域得到了广泛应用。 如需进一步了解傅里叶变换原理，可以参考相关资料。
　　快速傅里叶变换
　　计算机只能处理离散信号，使用离散傅里叶变换(DFT) 是计算机分析信号的基本方法。但是离散傅里叶变换的缺点是：计算量大，时间复杂度太高，当采样点数太高的时候，计算缓慢，由此出现了DFT的快速实现，即快速傅里叶变换FFT。
　　快速傅里叶变换（FFT）是计算量更小的离散傅里叶变换的一种实现方法，其逆变换被称为快速傅里叶逆变换（IFFT）。
　　示例
　　先对数据进行fft变换，然后再ifft逆变换。 import numpy as np #从fftpack中导入fft(快速傅里叶变化)和ifft(快速傅里叶逆变换)函数 from scipy.fftpack import fft,ifft  #创建一个随机值数组 x = np.array([1.0, 2.0, 1.0, -1.0, 1.5])  #对数组数据进行傅里叶变换 y = fft(x) print(＂fft: ＂) print(y) print(＂ ＂)  #快速傅里叶逆变换 yinv = ifft(y) print(＂ifft: ＂) print(yinv) print(＂ ＂)
　　输出 fft: [ 4.5       +0.j          2.08155948-1.65109876j -1.83155948+1.60822041j  -1.83155948-1.60822041j  2.08155948+1.65109876j]   ifft: [ 1. +0.j  2. +0.j  1. +0.j -1. +0.j  1.5+0.j]
　　可以看到fft，ifft返回的都是复数。ifft返回的结果中，复数的虚部都是0，实部与原始数据x一致。
　　这些点的频率无法计算，因为没有设置这N个点的时间长度。如不理解，不必深究，后面会介绍。 理解fft变换结果
　　我们知道，傅里叶变换把时域信号变为频域信号。在离散傅里叶变换中，频域信号由一系列不同频率的谐波（频率成倍数）组成。fft返回值是一个复数数组，每个复数表示一个正弦波。通常一个波形由振幅，相位，频率三个变量确定，可以从fft的返回值里，获取这些信息。
　　假设a是时域中的周期信号，采样频率为Fs，采样点数为N。如果A[N] = fft(a[N])，返回值A[N]是一个复数数组，其中： A[0]表示频率为0hz的信号，即直流分量。 A[1:N/2]包含正频率项，A[N/2:]包含负频率项。正频率项就是转化后的频域信号，通常我们只需要正频率项，即前面的n/2项，负频率项是计算的中间结果（正频率项的镜像值）。 每一项的频率计算：假设A[i]为数组中的元素，表示一个波形，该波形的频率 = i * Fs / N A[i] = real + j * imag，是一个复数，相位就是复数的辐角，相位 = arg(real/imag) 类似的，振幅就是复数的模，振幅 = sqrt(real^2+imag^2)。但是fft的返回值的模是放大值，直流分量的振幅放大了N倍，弦波分量的振幅放大了N/2倍。
　　频率分辨率
　　频率分辨率是离散傅里叶变换(DFT)频域相邻刻度之间的实际频率之差。采样时，数据采样了T秒(T = 采样点数N / 采样频率Fs)，信号的成分中周期最大也就是T秒，最低频率即＂基频＂就等于1 / T，也就是Fs / N，这就是频率分辨率。基频 = Fs / N，各个谐波的频率就是 i * Fs / N，这个公式用于计算各个波形的频率。
　　示例 import numpy as np from scipy.fftpack import fft  # 采样点数 N = 4000  # 采样频率 (根据采样定理，采样频率必须大于信号最高频率的2倍，信号才不会失真) Fs = 8000 x = np.linspace(0.0, N/Fs, N)  # 时域信号，包含：直流分量振幅1.0，正弦波分量频率100hz/振幅2.0, 正弦波分量频率150Hz/振幅0.5/相位np.pi y = 1.0 + 2.0 * np.sin(100.0 * 2.0*np.pi*x) + 0.5*np.sin(150.0 * 2.0*np.pi*x + np.pi)  # 进行fft变换 yf = fft(y)  # 获取振幅，取复数的绝对值，即复数的模 abs_yf = np.abs(yf)  # 获取相位，取复数的角度 angle_y=np.angle(yf)  # 直流信号 print(＂ 直流信号＂) print(＂振幅:＂, abs_yf[0]/N) # 直流分量的振幅放大了N倍  # 100hz信号 index_100hz = 100 * N // Fs # 波形的频率 = i * Fs / N，倒推计算索引：i = 波形频率 * N / Fs print(＂ 100hz波形＂) print(＂振幅:＂, abs_yf[index_100hz] * 2.0/N) # 弦波分量的振幅放大了N/2倍 print(＂相位:＂, angle_y[index_100hz])  # 150hz信号 index_150hz = 150 * N // Fs # 波形的频率 = i * Fs / N，倒推计算索引：i = 波形频率 * N / Fs print(＂ 150hz波形＂) print(＂振幅:＂, abs_yf[index_150hz] * 2.0/N) # 弦波分量的振幅放大了N/2倍 print(＂相位:＂, angle_y[index_150hz]) print(＂100hz与150hz相位差:＂,  angle_y[index_150hz] - angle_y[index_100hz]) print(＂ ＂)
　　输出 直流信号 振幅: 1.0  100hz波形 振幅: 1.9989359813189005 相位: -1.5315264186250062  150hz波形 振幅: 0.5008489983048182 相位: 1.6297011890497097 100hz与150hz相位差: 3.161227607674716
　　可以看到，正弦波的相位不一定从0开始，但波形之间的相位差确实s约等于一个pi(值跟采样频率与采样点数有关系)。 离散余弦变换(DCT)
　　由于许多要处理的信号都是实信号，在使用FFT时，对于实信号，傅立叶变换的共轭对称性导致在频域中有一半的数据冗余。
　　离散余弦变换（DCT）是对实信号定义的一种变换，变换后在频域中得到的也是一个实信号，相比离散傅里叶变换DFT而言, DCT可以减少一半以上的计算。DCT还有一个很重要的性质（能量集中特性）：大多书自然信号（声音、图像）的能量都集中在离散余弦变换后的低频部分，因而DCT在（声音、图像）数据压缩中得到了广泛的使用。由于DCT是从DFT推导出来的另一种变换，因此许多DFT的属性在DCT中仍然是保留下来的。
　　SciPy.fftpack中，提供了离散余弦变换(DCT)与离散余弦逆变换(IDCT)的实现。
　　示例 import numpy as np from scipy.fftpack import dct,idct y = dct(np.array([4., 3., 5., 10., 5., 3.])) print(y)
　　输出 [ 60.          -3.48476592 -13.85640646  11.3137085    6.   -6.31319305]
　　离散余弦逆变换(idct)，是离散余弦变换(DCT)的反变换。
　　示例 import numpy as np from scipy.fftpack import dct,idct y = idct(np.array([4., 3., 5., 10., 5., 3.])) print(y)
　　输出 [ 39.15085889 -20.14213562  -6.45392043   7.13341236   8.14213562   -3.83035081]