有趣的傅里叶变换与小波变换对比(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等。

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:/a/714170.html

如若内容造成侵权/违法违规/事实不符,请联系我们进行投诉反馈qq邮箱809451989@qq.com,一经查实,立即删除!

相关文章

Win11安装WSA 安卓系统,然后再电脑安装APK文件

参考文章&#xff1a; https://blog.csdn.net/m0_56076343/article/details/122334759 https://blog.csdn.net/u012514495/article/details/120885242 在微软的网站下载 打开&#xff1a;https://store.rg-adguard.net/ &#xff0c;如下图&#xff1a; 在 1 的那个地方&am…

二维数组与指针【C语言】

二维数组与指针 一维数组一维数组与指针二维数组二维数组与指针总结补充判断以下方式是否正确打印二维数组一维数组 int arr[] = {11, 22, 33, 44};arr:首地址(第一个元素的地址) 一维数组与指针 int arr[] = {11, 22, 33, 44};因为,arr表示的是首地址,等价于 int* p =…

谷粒商城实战(033 业务-秒杀功能4-高并发问题解决方案sentinel 2)

Java项目《谷粒商城》架构师级Java项目实战&#xff0c;对标阿里P6-P7&#xff0c;全网最强 总时长 104:45:00 共408P 此文章包含第332p-第p335的内容 熔断降级 开启对Feign远程服务的熔断保护机制 feign.sentinel.enabletrue 这里我们只是调用方加就行 被调用方不用加 正常…

PD19 Parallels Desktop 虚拟机 安装Windows10系统 操作步骤(保姆级教程,轻松上手)

Mac分享吧 文章目录 效果一、准备工作二、开始安装1、打开pd 19 虚拟机&#xff0c;点击右上角文件&#xff0c;新建2、通过下载好的镜像安装Windows10系统。找到镜像文件位置&#xff0c;安装&#xff0c;配置2、显示安装完成&#xff0c;打开Windows10系统 三、运行测试1、打…

SQLAlchemy:filter()和filter_by()的微妙差异

哈喽&#xff0c;大家好&#xff0c;我是木头左&#xff01; 在Python编程中&#xff0c;SQLAlchemy是一个强大的ORM&#xff08;对象关系映射&#xff09;工具&#xff0c;它允许使用Python代码来操作数据库。然而&#xff0c;对于新手来说&#xff0c;SQLAlchemy中的一些函数…

cocos开发的时候 wx.onShow在vscode里面显示红色

这个函数是在微信小游戏平台才会用到。 cocos识别不到wx这个变量。 可以改成下面的写法。 只要在变量前面加一个globalThis.就能识别这个变量了。也不报错了。 搞死强迫症了。orz 欢迎大家来玩我的微信小游戏。多多提意见啊。

【Java】图的初识

文章目录 【Java】图的初识图是什么图的基本组成部分图的类型图的表示方法图的常见操作 Java中图的表示方法邻接矩阵邻接表 常见操作图的遍历深度优先搜索&#xff08;DFS&#xff09;广度优先搜索&#xff08;BFS) 结论 【Java】图的初识 图是什么 图是一种数学概念&#xf…

[Linux] TCP协议介绍(1): TCP协议 数据格式、可靠性的控制、标记位... 简单介绍

上一篇文章, 针对UDP协议的格式、数据等内容做了一些简单的介绍. 并且提到, 在网络协议栈TCP/IP模型的传输层中, 有两个最具代表性的协议: UDP和TCP 下面就简单介绍分析一下TCP协议 TCP协议, 完整的称呼其实叫: 传输控制协议(Transmission Control Protocol) 从名字就可以看出…

Vue51-插件

一、插件的定义 vue里面的插件&#xff0c;类似于游戏的外挂。 vue中插件的本质&#xff1a;一个对象&#xff0c;里面必须包含install方法。 二、插件的使用 2-1、创建一个插件js文件&#xff08;写在src中plugins.js&#xff09; 2-2、应用插件&#xff1a;Vue.use(插件) …

基于Python+OpenCV高速公路行驶车辆的速度检测系统

简介&#xff1a; 基于Python和OpenCV的高速公路行驶车辆的速度检测系统旨在实时监测高速公路上的车辆&#xff0c;并测量它们的速度。该系统可以用于交通监控、道路安全管理等领域&#xff0c;为相关部门提供重要的数据支持。 系统实现&#xff1a; 视频流输入&#xff1a;系…

Python学习打卡:day07

day7 笔记来源于&#xff1a;黑马程序员python教程&#xff0c;8天python从入门到精通&#xff0c;学python看这套就够了 目录 day753、列表的常用操作课后练习题54、列表的循环遍历列表的遍历—— while 循环列表的遍历—— for 循环while 循环和 for 循环的对比练习 55、元组…

解决使用Jmeter进行测试时出现“302“,‘‘401“等用户未登录的问题

使用 JMeter 压力测试时解决登录问题的两种方法 在使用 JMeter 进行压力测试时&#xff0c;可能会遇程序存在安全验证&#xff0c;必须登录后才能对里面的具体方法进行测试&#xff1a; 如果遇到登录问题&#xff0c;通常是因为 JMeter 无法模拟用户的登录状态&#xff0c;导…

工程设计问题---行星轮系设计问题

该问题的主要目标是使汽车传动比的最大误差最小。为了使最大误差最小&#xff0c;对自动行星传动系统的齿轮齿数进行了计算。该问题包含6个整数变量和11个不同几何约束和装配约束的约束。 参考文献&#xff1a; Abhishek Kumar, Guohua Wu, Mostafa Z. Ali, Rammohan Malliped…

在 C++ 中使用不同平台的时间函数及比较

在 C 编程中&#xff0c;时间函数的选择对于性能测量、任务调度和时间戳记录至关重要。不同的操作系统提供了不同的时间函数&#xff0c;同时在同一个平台上&#xff0c;也可能有多种不同的时间函数可供选择。本文将介绍在 C 中常用的时间函数&#xff0c;并比较它们在不同平台…

【阿里云】2024年5月21日【算法岗暑期实习】面试经验分享

【阿里云】2024年5月21日【算法岗暑期实习】面试经验分享 面试流程&#xff1a;1个小时多。介绍WordEmbedding线性回归和逻辑回归的区别介绍 Information Gain&#xff08;信息增益&#xff09;IG做情感分析的时候存在多义词应该如何解决&#xff1f; 解题思路一&#xff1a;解…

【机器学习】机器学习重要分支——强化学习:从理论到实践

文章目录 强化学习&#xff1a;从理论到实践引言第一章 强化学习的基本概念1.1 什么是强化学习1.2 强化学习的基本组成部分1.3 马尔可夫决策过程 第二章 强化学习的核心算法2.1 Q学习2.2 深度Q网络&#xff08;DQN&#xff09;2.3 策略梯度方法 第三章 强化学习的应用实例3.1 游…

aarch64系统中给qt编译的可执行程序加上图标

在qt中编译的可执行程序图标是默认是一种格式 而且无法替换改图标&#xff0c;看起来不够美观&#xff0c;可以使用.desktop的链接文件来链接编译的执行程序&#xff0c;通过点击.desktop图标来去运行可执行程序。 实现步骤&#xff1a; 创建一个以.desktop结尾的文件并记得给…

ubuntu第三方库离线安装包(.deb离线安装方法;apt离线安装;离线安装deb)(docker离线安装、安装docker安装)

文章目录 方法1&#xff1a;Ubuntu Packages 网站下载离线包&#xff08;失败了&#xff0c;找不到包的可下载源&#xff0c;有的包有&#xff0c;有的包没有&#xff0c;不知道怎么回事&#xff09;操作步骤1. 在有网络的环境中&#xff0c;打开浏览器并访问 Ubuntu Packages …

【JS重点16】对象原型

目录 一&#xff1a;对象原型是什么 二&#xff1a;对象原型作用 三&#xff1a;constructor属性 四&#xff1a;如何赚钱 一&#xff1a;对象原型是什么 每个对象都有一个属性__proto__(称为原型对象),该属性是一个对象 __proto__是JS非标准属性在实例对象中&#xff0c;…

模型部署-什么是模型部署?

模型部署-什么是模型部署&#xff1f; 0总结1 什么是模型部署&#xff1f;1 部署流程&#xff1a;2 为什么模型部署这么复杂&#xff1f;3 模型部署学习建议&#xff1a; 2 详细介绍:1 模型转换2 模型优化3 模型压缩&#xff1a;剪枝&#xff1a;蒸馏稀疏化量化&#xff1a; 4 …