快速傅里叶变换（FFT）,离散傅里叶变换（DFT）
傅里叶变换
离散傅里叶变换（DFT）用于将一般时间序列变换到频域（即使它们是非周期的），即计算频谱。
计算DFT所需的DFT长度LDFT使用快速傅里叶变换（FFT）算法。np.fft.fft（x，L_DFT）必须选择得比信号x长。通常，如果DFT长度LDFT是2的幂，则fft算法最快。因此，在第12-14行中，检查DFT长度，发现L_DFT=2**14足够长，即比信号x（详见下文）长。
fs=8000 # sampling frequencyt=np.arange(0,2,1/fs) # time vectorf1=1230# fundamental frequency in Hzf2=1800# fundamental frequency in Hzsin1=np.sin(2*np.pi*f1*t)sin2=np.sin(2*np.pi*f2*t+np.pi)/2x=sin1+sin2# DFT of the signal mixtureL_DFT=2**14 # DFT length# we first check if the DFT length is longer than the signal lengthprint('Signal length is ' + str(len(x)))print('DFT length is ' + str(L_DFT))# perform DFT using the FFT algorithmX = np.fft.fft(x,L_DFT)# create a frequency vector to order negative /positive frequencies f1=np.linspace(0,fs/2,int(L_DFT/2)) # the positive frequencies (up to fs/2)f2=np.linspace(-fs/2,0,int(L_DFT/2)) # the negative frequenciesf = np.concatenate([f1, f2])# plot signalsplt.figure(figsize=(12,8))plt.subplot(2,2,1)plt.plot(t[:160],x[:160])plt.xlabel('time $t$ in seconds')plt.ylabel('$x(t)$')plt.subplot(2,2,2)plt.plot(f, np.abs(X))plt.title('Absolute value of the DFT of a signal mixture')plt.xlabel('frequency $f$ in Hz')plt.ylabel('$|x(f)|$')plt.subplot(2,2,3)plt.plot(f, np.real(X))plt.xlabel('frequency $f$ in Hz')plt.ylabel('$\mathrm{Re}\{x(f)\}$')plt.subplot(2,2,4)plt.plot(f, np.imag(X))plt.xlabel('frequency $f$ in Hz')plt.ylabel('$\mathrm{Im}\{x(f)\}$')plt.tight_layout()
sin2=np.sin(2np.pif2*t+np.pi)/2 why it add np.pi then divide 2 is that the way to discribe phi?np.real(X)np.imag(X)
可以看出，时域表示（左上角）不容易解释。频谱的大小（右上角）最清楚地表明，信号由两个正弦波部分组成，还可以看到信号在哪些频率下包含功率。从下方面板中的实部和虚部，我们可以看到光谱含量在实部和虚部之间“区分”。对于这个例子，我们可以得出结论，震级|x（f）|最容易解释。
提示：使用numpy函数fft可以简化上面示例中频率向量f的相当复杂的构造。fftshift（有关链接，请参见上面的ILO）。
找到合适的DFT/FFT长度-Helper函数计算下一次2的幂
DFT/FFT的一个重要参数是长度LDFT或LFFT。为了确定合适的FFT长度LFFT，我们感兴趣的是一个大于时域（输入）序列长度的数字。FFT通常被实现为各种类型和长度的输入信号有效计算DFT的众多算法的集合，通常使用软件库FFTW（这里为真正感兴趣的人撰写的科学论文）。在运行时，为信号类型（例如，实值、复值等）和长度LFFT选择有效算法。
2的幂有利于使用高效的Cooley–Tukey FFT算法或Radix-2实现FFT（）。因此，我们将使用下面的助手函数来确定信号长度的下一个二次幂。这个函数的确切工作方式超出了本课程的范围，但您应该能够使用它。
计算Next-Power-of-2的函数
以下函数nextPowerOf2（L）和nextPowerPOf2_simple（L）是计算适当的FFT长度L_FFT的两个实现，该长度是2的幂，并且总是大于信号长度L。如果您感兴趣，可以使用不同的L，例如，使用下面代码中的主注释中提供的示例。这些函数的确切工作方式超出了本实验表的范围。
def nextPowerOf2(L):'''Calculates the smallest power of 2 which is bigger than the length of input variable LThis helper function can be used to calculate an appropriate length for an DFT which is longer than the signal length L and is a power of 2.Input:L: intsignal lengthOutput:p: integer which is greater or equal than L and a power of 2Examples:for L in range(20):print('nextPowerOf2(L) for L='+str(L)+' is '+str(nextPowerOf2(L)))x=ones(11)L_FFT=nextPowerOf2(len(x))'''if (L<2):return 2# If n is a power of 2 then return n if (L and not(L & (L - 1))):return L# If n is NOT a power of 2 p = 1while (p < L) :p <<= 1 return pdef nextPowerOf2_simple(L):'''Calculates the smallest power of 2 which is bigger than the length of input variable LThis helper function can be used to calculate an appropriate length for an DFT which is longer than the signal length L and is a power of 2.Input:L: intsignal lengthOutput:p: integer which is greater or equal than n and a power of 2Examples:for L in range(20):print('nextPowerOf2(L) for L='+str(L)+' is '+str(nextPowerOf2_simple(L)))x=ones(11)L_FFT=nextPowerOf2_simple(len(x))'''return int(np.max([2,2**np.ceil(np.log2(L))]))for L in range(20):print('nextPowerOf2(L) for L='+str(L)+' is '+str(nextPowerOf2(L)))# this version would create a warning for L=0, to put in 0 as length of the signal should also not happenfor L in range(20):print('nextPowerOf2_simple(L) for L='+str(L)+' is '+str(nextPowerOf2_simple(L)))
