有趣的傅里叶变换与小波变换对比（Python）

不严谨的说，时域和频域分析就是在不同的空间看待问题的，不同空间所对应的原子(基函数)是不同的。你想一下时域空间的基函数是什么？频域空间的基函数是什么？一般的时-频联合域空间的基函数是什么？小波域空间的基函数是什么？

有的空间域比较容易分析，有的空间域不容易分析。

举个例子吧，首先加载一个双曲Chirp信号，数据的采样频率为2048Hz，第一个Chirp信号持续时间为0.1~0.68秒，第二个Chirp信号持续时间为0.1~0.75 秒，第一个Chirp信号在时间t处的瞬时频率为（单位Hz）：

第二个Chirp信号在时间t处的瞬时频率为（单位Hz）：

看一下从时域空间看待的时域图

然后看一下频域空间的频谱图
傅里叶变换（FT）比较擅长识别信号中存在的频率分量， 但是FT无法定位频率分量。绘制上面信号的幅值谱，并放大0到200Hz之间的区域

再看一下一般的时频域空间的时频谱图，以短时傅里叶变换为例
傅里叶变换不提供时间信息，为了定位频率，短时傅里叶变换STFT方法将信号分割成不同的窗，并对每个窗执行FT。STFT的时频分析窗口如下：

STFT提供了信号时间-频率域中的一些信息， 但是选择窗的大小是关键。对于STFT时频分析，选择更短的窗以牺牲频率分辨率为代价从而获得良好的时间分辨率。相反，选择较大的窗以时间分辨率为代价从而获得良好的频率分辨率（著名的测不准原理）。一旦STFT的分析窗确定后，将在整个分析中保持不变（最致命的缺陷）。以 200 毫秒的时间窗大小绘制上述双曲Chirp信号的频谱图，频谱图上的瞬时频率为黑色虚线段。

然后绘制时间窗大小为50毫秒的频谱图

两个图的结果是显而易见的，没有单一的窗口大小可以解析此类信号的整个频率信息。
最后看一下小波空间对应的小波时频谱图
连续小波变换 CWT是为了克服 STFT中固有的时频分辨率问题。CWT的时频分析窗口如下：

CWT和人类的听觉系统非常一致:在低频处有更好的频率定位能力，在高频处有更好的时间定位能力。绘制 CWT时尺度谱（尺度谱是作为时间和频率绘制的 CWT的绝对值），因为CWT 中的频率是对数的，所以使用对数频率轴。

从图中可以清楚地看出信号中两个双曲Chirp信号的存在，CWT可以比较准确估计持续时间的瞬时频率，而无需担心选择窗的大小。要了解小波系数幅度增长速度有多快，可以看一下3-D 图

在尺度谱上绘制一下瞬时频率，可见瞬时频率与尺度谱特征非常吻合

看到了吧，从不同的空间域（角度）看待问题，分析的难度也不一样。

开始正题Wavelet vs Fourier transform

#Import necessary libraries
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns


import pywt
from scipy.ndimage import gaussian_filter1d
from scipy.signal import chirp
import matplotlib.gridspec as gridspec
from scipy import signal
from skimage import filters,img_as_float
from skimage.io import imread, imshow
from skimage.color import rgb2hsv, rgb2gray, rgb2yuv
from skimage import color, exposure, transform
from skimage.exposure import equalize_hist
from scipy import fftpack, ndimage
t_min=0
t_max=10
fs=100
dt = 1/fs
time = np.linspace(t_min, t_max, 1500)
#To understand the behaviour of scale, we used a smooth constant signal with a discontinuity. Adding discontinuity to the constant will have a rectangular shape.
w = chirp(time, f0=10, f1=50, t1=10, method='quadratic')


#Compute Wavelet Transform
scale = [10,20,30,50,100]


#Plot signal, FFT, and scalogram(to represent wavelet transform)
fig,axes =  plt.subplots(nrows=1,ncols=5,figsize=(25,4))
for i in range(2):
  for j in range(5):
    #Scalogram
    scales = np.arange(1,scale[j],1)
    coef,freqs = pywt.cwt(w,scales,'morl')
    freqs = pywt.scale2frequency('morl',scales,precision=8)
    if i == 0:
      axes[j].set_title("Scalogram from scale {} to {}".format(1,scale[j]))
    if i == 0:
      axes[j].pcolormesh(time, scales, coef,cmap='Greys')
      axes[j].set_ylabel("Scale")
plt.show();

scales = np.arange(1,20,1)
coef,freqs = pywt.cwt(w,scales,'morl',1/fs)
fig,axes =  plt.subplots(nrows=1,ncols=2,figsize=(12,5))
axes[0].set_title("Scalogram")
axes[0].pcolormesh(time, scales, coef,cmap='Greys')
axes[0].set_xlabel("Time")
axes[0].set_ylabel("Scale")
axes[1].set_title("Spectrogram")
axes[1].pcolormesh(time, freqs, coef,cmap='Greys')
axes[1].set_xlabel("Time")
axes[1].set_ylabel("Pseudo Frequency")
plt.show();

families = ['gaus1','gaus2','gaus3','gaus4','gaus5','gaus6','gaus7','gaus8','mexh','morl']
cols = 5
rows = 4
scales = np.arange(1,20,1)
fig,axes =  plt.subplots(nrows = rows,ncols=5,figsize=(3*cols,2*rows))
fig.tight_layout(pad=1.0, w_pad=1.0, h_pad=3)
for i,family in enumerate(families):
  c = i%5
  r = round(i//5)
  coef,freqs = pywt.cwt(w,scales,family,1/fs)
  psi, x = pywt.ContinuousWavelet(family).wavefun(level=10)
  axes[r*2,c].set_title(family)
  axes[(r*2)+1,c].pcolormesh(time, freqs, coef,cmap='Blues')
  axes[(r*2)+1,c].set_xlabel("Time")
  axes[(r*2)+1,c].set_ylabel("Scale")
  axes[r*2,c].plot(x, psi)
  axes[r*2,c].set_xlabel("X")
  axes[r*2,c].set_ylabel("Psi")

fs = 100 #Sampling frequency
time = np.arange(-3,3,1/fs) #create time
n = len(time)
T=1/fs
print("We consider {} samples".format(n))
constant = np.ones(n) #Amblitude will be one(constant value)
freq =  np.linspace(-1.0/(2.0*T), 1.0/(2.0*T), n)


#Compute Fourier transform of Constant signal
fft = fftpack.fft(constant)
freq = fftpack.fftfreq(time.shape[0],T)
phase  = np.angle(fft)
phase  = phase / np.pi


#Compute Wavelet Transform
scales = np.arange(1,6,1)
coef,freqs = pywt.cwt(constant,scales,'gaus1')


#Plot signal, FFT, and scalogram(to represent wavelet transform)
fig,axes =  plt.subplots(ncols=3,figsize=(18,4))


#Signal
axes[0].set_title("Constant")
axes[0].plot(time, constant)
axes[0].set_xlabel("Time")
axes[0].set_ylabel("Amplitude")


#Fourier
axes[1].set_title("Fourier Transform")
axes[1].plot(freq, np.abs(fft)/n)
axes[1].set_xlabel("Frequency")
axes[1].set_ylabel("Magnitude")


#Scalogram
axes[2].set_title("Scalogram")
axes[2].pcolormesh(time, scales, coef,cmap='bone')
axes[2].set_xlabel("Time")
axes[2].set_ylabel("Scale")
plt.show();

constant[300:340]=0


#Compute Fourier transform of Constant signal
fft = fftpack.fft(constant)
phase  = np.angle(fft)
phase  = phase / np.pi


#Compute Wavelet Transform
scales = np.arange(1,6,1)
coef,freqs = pywt.cwt(constant,scales,'gaus1')


#Plot signal, FFT, and scalogram(to represent wavelet transform)
fig,axes =  plt.subplots(ncols=3,figsize=(18,4))


#Signal
axes[0].set_title("Constant")
axes[0].plot(time, constant)
axes[0].set_xlabel("Time")
axes[0].set_ylabel("Amplitude")


#Fourier
axes[1].set_title("Fourier Transform")
axes[1].plot(freq, np.abs(fft)/n)
axes[1].set_xlabel("Frequency")
axes[1].set_ylabel("Magnitude")


#Scalogram
axes[2].set_title("Scalogram")
axes[2].pcolormesh(time, scales, coef,cmap='bone')
axes[2].set_xlabel("Time")
axes[2].set_ylabel("Scale")
plt.show();

N = 50000 #number of samples
fs = 1000 #sample frequency
T = 1/fs #interval
time = np.linspace(-(N*T), N*T, N)
rect = np.zeros(time.shape)
for i in range(time.shape[0]):
    if time[i] > -0.5 and time[i] < 0.5:
        rect[i] = 1.0
print("We consider {} samples".format(N))
freq =  np.linspace(-1.0/(2.0*T), 1.0/(2.0*T), N)


#compute Fourier Trainsform
fft_rect = np.fft.fft(rect)
fr = np.fft.fftfreq(N)
phase  = np.angle(fft_rect)
phase  = phase / np.pi
freqrect = np.fft.fftfreq(time.shape[-1])
fft_rect = np.fft.fftshift(fft_rect)


#compute wavelet transform
scales = np.arange(1,25,1)
coef,freqs = pywt.cwt(rect,scales,'gaus1')


#Plot
#signal
fig,axes =  plt.subplots(ncols=3,figsize=(21,5))
axes[0].set_title("Rectangular signal")
axes[0].plot(time, rect)
axes[0].set_xlim(-1,1)
axes[0].set_xlabel("Time")
axes[0].set_ylabel("rectangular pulse")


#Fourier transform
axes[1].set_title("Fourier Transform")
axes[1].plot(freq,np.abs(fft_rect)*2/fs)
axes[1].set_xlim(-40,40)
axes[1].set_xlabel("Frequency")
axes[1].set_ylabel("Magnitude")


#wavelet
axes[2].set_title("Scalogram ")
axes[2].pcolormesh(time, scales, coef,cmap='bone')
axes[2].set_xlim(-2,2)
axes[2].set_xlabel("Time")
axes[2].set_ylabel("Scale")
plt.show();

fs = 1000 #sampling frequency
interval = 1/fs #sampling interval
t_min = -1 #start time
t_max = 1 # end time
dt=1/fs
time = np.arange(t_min,t_max,interval)


n = len(time)
print("We consider {} samples".format(n))


f = (fs/2)*np.linspace(0,1,int(n/2)) #frequency


freq = [200,130] #signal frequencies
scales1 = np.arange(1,20,1)


#Create signal with 200 hz frequency
sinewave1 = np.sin(2*np.pi*freq[0]*time)
new = sinewave1/np.square(time)
#compute fourier transform
fft1 = np.fft.fft(sinewave1)
fr = np.fft.fftfreq(n, d=dt)
phase  = np.angle(fft1)
phase  = phase / np.pi
fft1 = fft1[0:int(n/2)]


#compute wavelet
coef1,freqs1 = pywt.cwt(sinewave1,scales1,'morl')


#plot
gs = gridspec.GridSpec(2,2)
gs.update(left=0, right=4,top=2,bottom=0, hspace=.2,wspace=.1)
ax = plt.subplot(gs[0, :])
ax.set_title("Sinusoidal Signal - 200 Hz")
ax.plot(time,sinewave1)
ax.set_xlabel("Time(s)")
ax.set_ylabel("Amplitude")
ax2 = plt.subplot(gs[1, 0])
ax2.plot(f,np.abs(fft1)*2/fs)
ax2.set_xlabel("Frequency")
ax2.set_ylabel("DFT values")
ax3 = plt.subplot(gs[1, 1])
ax3.pcolormesh(time, freqs1/dt, coef1)
ax3.set_xlabel("Time")
ax3.set_ylabel("Frequency")
plt.show;

scales = np.arange(1,20,1)
#Create signal with 130 hz frequency
sinewave2 = np.sin(2*np.pi*freq[1]*time)


#compute fourier transform
fft2 = np.fft.fft(sinewave2)
fft2=fft2[0:int(n/2)]


#compute wavelet
coef2,freqs2 = pywt.cwt(sinewave2,scales,'morl')
#plot
gs = gridspec.GridSpec(2,2)
gs.update(left=0, right=4,top=2,bottom=0, hspace=.2,wspace=.1)
ax = plt.subplot(gs[0, :])
ax.set_title("Sinusoidal Signal - 130 Hz")
ax.plot(time,sinewave2)
ax.set_xlabel("Time(s)")
ax.set_ylabel("Amplitude")
ax2 = plt.subplot(gs[1, 0])
ax2.plot(f,np.abs(fft2)*2/fs)
ax2.set_xlim(0,300)
ax2.set_xlabel("Frequency")
ax2.set_ylabel("Magnitude")
ax3 = plt.subplot(gs[1, 1])
ax3.pcolormesh(time, freqs2/dt, coef2)
ax3.set_xlabel("Time")
ax3.set_ylabel("Scale")
plt.show();

scales = np.arange(1,30,1)


sum = sinewave1+sinewave2
fft3 = np.fft.fft(sum)
fft3=fft3[0:int(n/2)]
coef3,freqs3 = pywt.cwt(sum,scales,'morl')


#plot
gs = gridspec.GridSpec(2,2)
gs.update(left=0, right=4,top=2,bottom=0, hspace=.2,wspace=.1)




ax = plt.subplot(gs[0, :])
ax.set_title("Sum of Sinusoidal")
ax.plot(time,sum)
ax.set_xlabel("Time(s)")
ax.set_ylabel("Amplitude")
ax2 = plt.subplot(gs[1, 0])
ax2.plot(f,np.abs(fft3)*2/fs)
ax2.set_xlabel("Frequency")
ax2.set_ylabel("Magnitude")
ax3 = plt.subplot(gs[1, 1])
ax3.pcolormesh(time, freqs3/dt, coef3)
ax3.set_xlabel("Time")
ax3.set_ylabel("Frequency")
plt.show();

scales = np.arange(1,40,1)
#creating two cosine waves at 50 hz frequency
coswave = np.cos(2*np.pi*50*time)
#compute Fourier Trainsform of amblitude with 10 Hz
fft = np.fft.fft(coswave)
fft=fft[0:int(n/2)]
coef,freqs = pywt.cwt(coswave,scales,'morl')
#plot
gs = gridspec.GridSpec(2,2)
gs.update(left=0, right=4,top=2,bottom=0, hspace=.2,wspace=.1)




ax = plt.subplot(gs[0, :])
ax.set_title("Cosine Signal - 50 Hz")
ax.plot(time,coswave)
ax.set_xlabel("Time(s)")
ax.set_ylabel("Amplitude")
ax2 = plt.subplot(gs[1, 0])
ax2.plot(f,np.abs(fft)*2/fs)
ax2.set_xlabel("Frequency")
ax2.set_ylabel("Magnitude")
ax3 = plt.subplot(gs[1, 1])
ax3.pcolormesh(time, freqs/dt, coef)
ax3.set_xlabel("Time")
ax3.set_ylabel("Frequency")
plt.show()

Inject a Signal to the Sum signal

scales = np.arange(1,40,1)
sum_disc = sum+np.real(np.exp(-50*(time-0.4)**2)*np.exp(1j*2*np.pi*400*(time-0.4)))
fft = np.fft.fft(sum_disc)
fft=fft[0:int(n/2)]
coef,freqs = pywt.cwt(sum_disc,scales,'morl')




#plot
gs = gridspec.GridSpec(2,2)
gs.update(left=0, right=4,top=2,bottom=0, hspace=.2,wspace=.1)




ax = plt.subplot(gs[0, :])
ax.set_title("Signal with variation")
ax.plot(time,sum_disc)
ax.set_xlabel("Time(s)")
ax.set_ylabel("Amplitude")
ax2 = plt.subplot(gs[1, 0])
ax2.plot(f,np.abs(fft)*2/fs)
ax2.set_xlim(0,500)
ax2.set_xlabel("Frequency")
ax2.set_ylabel("Magnitude")
ax3 = plt.subplot(gs[1, 1])
ax3.pcolormesh(time, freqs*1000, coef)
ax3.set_xlabel("Time")
ax3.set_ylabel("Frequency")

Non stationary signals

size = len(time)//3
scales = np.arange(1,31,1)
sig = np.zeros(time.shape)
sig[:size]=np.sin(2*np.pi*200*time[:size])
sig[size:size*2]=np.sin(2*np.pi*130*time[size:size*2])
sig[size*2:]=np.cos(2*np.pi*50*time[size*2:])
fft = np.fft.fft(sig)
fft=fft[0:int(n/2)]
coef,freqs = pywt.cwt(sig,scales,'gaus8')
stft_f, stft_t, Sxx = signal.spectrogram(sig, fs,window='hann', nperseg=64)
#plot
gs = gridspec.GridSpec(2,2)
gs.update(left=0, right=4,top=2,bottom=0, hspace=.2,wspace=.1)




ax = plt.subplot(gs[0, 0])
ax.set_title("Sinusoidal Signal- Frequency vary over time")
ax.plot(time,sig)
ax.set_xlabel("Time(s)")
ax.set_ylabel("Amplitude")
ax1 = plt.subplot(gs[0, 1])
ax1.plot(f,np.abs(fft)*2/fs)
ax1.set_xlabel("Frequency")
ax1.set_ylabel("Magnitude")
ax2 = plt.subplot(gs[1, 0])
ax2.pcolormesh(stft_t, stft_f, Sxx)
ax2.set_ylabel("Frequency")
ax3 = plt.subplot(gs[1, 1])
ax3.pcolormesh(time, freqs/dt, coef)
ax3.set_xlabel("Time")
ax3.set_ylabel("Frequency")

Linear Chirp Signal

def plot_chirp_transforms(type_,f0,f1):
  #Create linear chirp signal with frequency between 50Hz and 10Hz
  t_min=0
  t_max=10
  time = np.linspace(t_min, t_max, 1500)
  N = len(time)
  interval = (t_min+t_max)/N
  fs = int(1/interval)
  dt=1/fs
  f = (fs/2)*np.linspace(0,1,int(N/2))
  w1 = chirp(time, f0=f0, f1=f1, t1=10, method=type_.lower())
  #Compute FFT
  w1fft = np.fft.fft(w1)
  w1fft=w1fft[0:int(N/2)]
  #Compute Wavelet transform
  scales=np.arange(1,50,1)
  wcoef,wfreqs = pywt.cwt(w1,scales,'morl')


  #Compute Short Time Fourier transfomr
  stft_f, stft_t, Sxx = signal.spectrogram(w1, fs,window='hann', nperseg=64,noverlap=32)


  #Plot the results
  gs = gridspec.GridSpec(2,2)
  gs.update(left=0, right=4,top=2,bottom=0, hspace=.2,wspace=.1)
  ax = plt.subplot(gs[0, 0])
  ax.set_title("Chirp - "+type_+" ({}Hz to {}Hz)".format(f0,f1))
  ax.plot(time,w1)
  ax.set_xlabel("Time(s)")
  ax.set_ylabel("Amplitude")
  ax1 = plt.subplot(gs[0, 1])
  ax1.plot(f,np.abs(w1fft)*2/fs)
  plt.grid()
  ax1.set_title("FFT - "+type_+" chirp signal ({}Hz to {}Hz)".format(f0,f1))
  ax1.set_xlabel("Frequency")
  ax1.set_ylabel("Magnitude")
  ax2 = plt.subplot(gs[1, 0])
  ax2.set_title("STFT - "+type_+" chirp signal ({}Hz to {}Hz)".format(f0,f1))
  ax2.pcolor(stft_t, stft_f, Sxx,cmap='copper')
  ax2.set_xlabel("Time")
  ax2.set_ylabel("Frequency")
  ax3 = plt.subplot(gs[1, 1])
  ax3.set_title("WT - "+type_+" chirp signal ({}Hz to {}Hz)".format(f0,f1))
  ax3.pcolor(time, wfreqs/dt, wcoef,cmap='copper')
  ax3.set_ylim(5,75)
  ax3.set_xlabel("Time")
  ax3.set_ylabel("Frequency")
plot_chirp_transforms('Linear',50,10)

plot_chirp_transforms('Linear',10,50)

plot_chirp_transforms('Logarithmic',50,10)

plot_chirp_transforms('Hyperbolic',50,10)

Triangular Signal

N = 1000 #number of samples
fs = 1000 #sample frequency
T = 1/fs #interval
time = np.linspace(-2, 2, N)
tri = np.where(np.abs(time)<=.5,1,0)
tri = np.where(tri==1,.5-np.abs(time),0)


print("We consider {} samples".format(N))
scales = np.arange(1,51,1)
coef,freqs = pywt.cwt(tri,scales,'gaus1')
#compute Fourier Trainsform
fft = np.fft.fft(tri)
freq = np.fft.fftfreq(time.shape[-1],T)
fftShift = np.fft.fftshift(fft)
freqShift=np.fft.fftshift(freq)






#Plot signal and FFT
fig,axes =  plt.subplots(nrows=2,ncols=3,figsize=(24,10))
axes[0,0].set_title("Triangular signal")
axes[0,0].plot(time, tri)
axes[0,0].set_xlabel("Time")
axes[0,0].set_ylabel("Triangular pulse")
axes[0,1].set_title("Fourier Transform")
axes[0,1].plot(freqShift,np.abs(fftShift)*2/fs)
axes[0,1].set_xlabel("Frequency")
axes[0,1].set_xlim(-50,50)
axes[0,1].set_ylabel("Magnitude")
axes[0,2].set_title("Scalogram")
axes[0,2].pcolor(time,scales,coef,cmap='copper')
axes[0,2].set_xlim(-2,2)
axes[0,2].set_xlabel("Time")
axes[0,2].set_ylabel("Scale")


indices = np.where(tri>0)[0]
new_indices = indices+120
temp = tri.copy()
tri = np.zeros(tri.shape)
tri[new_indices] = temp[indices]




coef,freqs = pywt.cwt(tri,scales,'gaus1')
#compute Fourier Trainsform
fft = np.fft.fft(tri)
freq = np.fft.fftfreq(time.shape[-1],T)
fftShift = np.fft.fftshift(fft)
freqShift=np.fft.fftshift(freq)


#Plot signal and FFT
axes[1,0].plot(time, tri)
axes[1,0].set_xlabel("Time")
axes[1,0].set_ylabel("Triangular pulse")
axes[1,1].plot(freqShift,np.abs(fftShift)*2/fs)
axes[1,1].set_xlim(-50,50)
axes[1,1].set_xlabel("Frequency")
axes[1,1].set_ylabel("Magnitude")
axes[1,2].pcolor(time,scales,coef,cmap='copper')
axes[1,2].set_xlim(-2,2)
axes[1,2].set_xlabel("Time")
axes[1,2].set_ylabel("Scale")

Audio

from scipy.io import wavfile
import scipy
#Read Audio and compute time
sr, data = wavfile.read('5.wav')
dt = 1/sr
time = np.arange(0,1,dt)
#Find FFT and frequencies
fft_aud = np.fft.fft(data)
fft_aud=fft_aud[0:int(sr/2)]
freq = (sr/2)*np.linspace(0,1,int(sr/2))
plt.plot(time,data)


#Compute STFT
stft_f, stft_t, Sxx = signal.spectrogram(data, sr,window='hann', nperseg=256,noverlap=64)


#Compute Wavelet Transform(morlet)
widths = np.arange(1, 31)
wt,wfreqs = pywt.cwt(data,widths,'morl')
gs = gridspec.GridSpec(2,2)
gs.update(left=0, right=4,top=2,bottom=0, hspace=.2,wspace=.1)


ax0 = plt.subplot(gs[0, 0])
ax0.set_title("Audio Signal")
ax0.plot(time,data)
ax0.set_xlabel("Time")
ax0.set_ylabel("Amplitude")
ax1 = plt.subplot(gs[0, 1])
ax1.plot(freq, np.abs(fft_aud))
ax1.set_xlabel("Frequency")
ax1.set_ylabel("Magnitude")
ax2 = plt.subplot(gs[1, 0])
ax2.pcolor(stft_t, stft_f, Sxx, cmap='copper')
ax2.set_xlabel("Time")
ax2.set_ylabel("Frequency")
ax3 = plt.subplot(gs[1, 1])
ax3.pcolormesh(time, wfreqs/dt, wt)
ax3.set_xlabel("Time")
ax3.set_ylabel("Frequency")
plt.show()

from PIL import Image
# open the original image
original_img = Image.open("parrot1.jpg")




#rotate image
rot_180 = original_img.rotate(180, Image.NEAREST, expand = 1)


# close all our files object


I = np.array(original_img)
I_rot = np.array(rot_180)




original_img.close()


I_grey = rgb2gray(I)
I_rot_grey = rgb2gray(I_rot)




fft2 = fftpack.fft2(I_grey)
fftshift = fftpack.fftshift(fft2)
fftrot2 = fftpack.fft2(I_rot_grey)
fftrotshift = fftpack.fftshift(fftrot2)


coeffs2 = pywt.dwt2(I_grey, 'haar')
cA, (cH, cV, cD) = coeffs2
titles = ['Approximation', ' Horizontal detail','Vertical detail', 'Diagonal detail']
coeffs3 = pywt.dwt2(I_rot_grey, 'haar')
cA1, (cH1, cV1, cD1) = coeffs3
fig,axes = plt.subplots(ncols=6,nrows=2,figsize=(24,8))


axes[0,0].set_title("Image")
axes[0,0].imshow(img_as_float(I_grey),cmap='gray')
axes[0,1].set_title("FFT")
axes[0,1].imshow(np.log(np.abs(fftshift)),cmap='gray')
axes[1,0].set_title("Flip Image")
axes[1,0].imshow(img_as_float(I_rot_grey),cmap='gray')
axes[1,1].set_title("FFT-flip image")
axes[1,1].imshow(np.log(np.abs(fftrotshift)),cmap='gray')




for idx,coef in enumerate((cA,cH,cV,cD)):
  axes[0,idx+2].set_title(titles[idx])
  axes[0,idx+2].imshow(coef,cmap='gray')
for idx,coef in enumerate((cA1,cH1,cV1,cD1)):
  axes[1,idx+2].set_title(titles[idx])
  axes[1,idx+2].imshow(coef,cmap='gray')
plt.show();

Trapezoid

N = 5000 #number of samples
fs = 1000 #sample frequency
T = 1/fs #interval
time = np.linspace(-5, 5, N)
trapzoid_signal = (time*np.where(time>0,1,0))-((time-1)*np.where((time-1)>0,1,0))-((time-2)*np.where((time-2)>0,1,0))+((time-3)*np.where((time-3)>0,1,0))


#tra = trapzoid_signal(time)


scales = np.arange(1,51,1)
coef,freqs = pywt.cwt(trapzoid_signal,scales,'gaus1')
#compute Fourier Trainsform
fft = np.fft.fft(trapzoid_signal)
freq = np.fft.fftfreq(time.shape[-1],T)
fftShift = np.fft.fftshift(fft)
freqShift=np.fft.fftshift(freq)


#Plot signal and FFT
fig,axes =  plt.subplots(nrows=2,ncols=3,figsize=(24,10))
axes[0,0].set_title("Trapezoidal signal")
axes[0,0].plot(time, trapzoid_signal)
axes[0,0].set_xlabel("Time")
axes[0,0].set_ylabel("Trapezoidal pulse")
axes[0,1].set_title("Fourier Transform - trapezoidal")
axes[0,1].plot(freqShift,np.abs(fftShift)*2/fs)
axes[0,1].set_xlim(-20,20)
axes[0,1].set_xlabel("Frequency")
axes[0,1].set_ylabel("Magnitude")
axes[0,2].set_title("Scalogram - trapezoidal")
axes[0,2].pcolor(time,scales,coef,cmap='BrBG')
axes[0,2].set_xlim(-5,5)
axes[0,2].set_xlabel("Time")
axes[0,2].set_ylabel("Scale")


trapzoid_signal = ((time+1)*np.where((time+1)>0,1,0))-(time*np.where(time>0,1,0))-((time-1)*np.where((time-1)>0,1,0))+((time-2)*np.where((time-2)>0,1,0))
coef,freqs = pywt.cwt(trapzoid_signal,scales,'gaus1')
#compute Fourier Trainsform
fft = np.fft.fft(trapzoid_signal)
freq = np.fft.fftfreq(time.shape[-1],T)
fftShift = np.fft.fftshift(fft)
freqShift=np.fft.fftshift(freq)


#Plot signal and FFT
axes[1,0].plot(time, trapzoid_signal)
axes[1,0].set_xlabel("Time")
axes[1,0].set_ylabel("Trapezoidal pulse")
axes[1,1].plot(freqShift,np.abs(fftShift)*2/fs)
axes[1,1].set_xlim(-20,20)
axes[1,1].set_xlabel("Frequency")
axes[1,1].set_ylabel("Magnitude")
axes[1,2].pcolor(time,scales,coef,cmap='BrBG')
axes[1,2].set_xlim(-5,5)
axes[1,2].set_xlabel("Time")
知乎学术咨询：
https://www.zhihu.com/consult/people/792359672131756032?isMe=1
axes[1,2].set_ylabel("Scale")

工学博士，担任《Mechanical System and Signal Processing》《中国电机工程学报》《控制与决策》等期刊审稿专家，擅长领域：现代信号处理，机器学习，深度学习，数字孪生，时间序列分析，设备缺陷检测、设备异常检测、设备智能故障诊断与健康管理PHM等。