有趣的傅里叶变换与小波变换对比(Python)

2024-06-16 14:52

本文主要是介绍有趣的傅里叶变换与小波变换对比(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 snsimport 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):#Scalogramscales = 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%5r = 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)) #frequencyfreq = [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 10Hzt_min=0t_max=10time = np.linspace(t_min, t_max, 1500)N = len(time)interval = (t_min+t_max)/Nfs = int(1/interval)dt=1/fsf = (fs/2)*np.linspace(0,1,int(N/2))w1 = chirp(time, f0=f0, f1=f1, t1=10, method=type_.lower())#Compute FFTw1fft = np.fft.fft(w1)w1fft=w1fft[0:int(N/2)]#Compute Wavelet transformscales=np.arange(1,50,1)wcoef,wfreqs = pywt.cwt(w1,scales,'morl')#Compute Short Time Fourier transfomrstft_f, stft_t, Sxx = signal.spectrogram(w1, fs,window='hann', nperseg=64,noverlap=32)#Plot the resultsgs = 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 objectI = 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等。

这篇关于有趣的傅里叶变换与小波变换对比(Python)的文章就介绍到这儿,希望我们推荐的文章对编程师们有所帮助!



http://www.chinasem.cn/article/1066744

相关文章

如何使用 Python 读取 Excel 数据

《如何使用Python读取Excel数据》:本文主要介绍使用Python读取Excel数据的详细教程,通过pandas和openpyxl,你可以轻松读取Excel文件,并进行各种数据处理操... 目录使用 python 读取 Excel 数据的详细教程1. 安装必要的依赖2. 读取 Excel 文件3. 读

Python的time模块一些常用功能(各种与时间相关的函数)

《Python的time模块一些常用功能(各种与时间相关的函数)》Python的time模块提供了各种与时间相关的函数,包括获取当前时间、处理时间间隔、执行时间测量等,:本文主要介绍Python的... 目录1. 获取当前时间2. 时间格式化3. 延时执行4. 时间戳运算5. 计算代码执行时间6. 转换为指

利用Python调试串口的示例代码

《利用Python调试串口的示例代码》在嵌入式开发、物联网设备调试过程中,串口通信是最基础的调试手段本文将带你用Python+ttkbootstrap打造一款高颜值、多功能的串口调试助手,需要的可以了... 目录概述:为什么需要专业的串口调试工具项目架构设计1.1 技术栈选型1.2 关键类说明1.3 线程模

Python ZIP文件操作技巧详解

《PythonZIP文件操作技巧详解》在数据处理和系统开发中,ZIP文件操作是开发者必须掌握的核心技能,Python标准库提供的zipfile模块以简洁的API和跨平台特性,成为处理ZIP文件的首选... 目录一、ZIP文件操作基础三板斧1.1 创建压缩包1.2 解压操作1.3 文件遍历与信息获取二、进阶技

Python Transformers库(NLP处理库)案例代码讲解

《PythonTransformers库(NLP处理库)案例代码讲解》本文介绍transformers库的全面讲解,包含基础知识、高级用法、案例代码及学习路径,内容经过组织,适合不同阶段的学习者,对... 目录一、基础知识1. Transformers 库简介2. 安装与环境配置3. 快速上手示例二、核心模

Python正则表达式语法及re模块中的常用函数详解

《Python正则表达式语法及re模块中的常用函数详解》这篇文章主要给大家介绍了关于Python正则表达式语法及re模块中常用函数的相关资料,正则表达式是一种强大的字符串处理工具,可以用于匹配、切分、... 目录概念、作用和步骤语法re模块中的常用函数总结 概念、作用和步骤概念: 本身也是一个字符串,其中

Python使用getopt处理命令行参数示例解析(最佳实践)

《Python使用getopt处理命令行参数示例解析(最佳实践)》getopt模块是Python标准库中一个简单但强大的命令行参数处理工具,它特别适合那些需要快速实现基本命令行参数解析的场景,或者需要... 目录为什么需要处理命令行参数?getopt模块基础实际应用示例与其他参数处理方式的比较常见问http

python实现svg图片转换为png和gif

《python实现svg图片转换为png和gif》这篇文章主要为大家详细介绍了python如何实现将svg图片格式转换为png和gif,文中的示例代码讲解详细,感兴趣的小伙伴可以跟随小编一起学习一下... 目录python实现svg图片转换为png和gifpython实现图片格式之间的相互转换延展:基于Py

Python中的getopt模块用法小结

《Python中的getopt模块用法小结》getopt.getopt()函数是Python中用于解析命令行参数的标准库函数,该函数可以从命令行中提取选项和参数,并对它们进行处理,本文详细介绍了Pyt... 目录getopt模块介绍getopt.getopt函数的介绍getopt模块的常用用法getopt模

Python利用ElementTree实现快速解析XML文件

《Python利用ElementTree实现快速解析XML文件》ElementTree是Python标准库的一部分,而且是Python标准库中用于解析和操作XML数据的模块,下面小编就来和大家详细讲讲... 目录一、XML文件解析到底有多重要二、ElementTree快速入门1. 加载XML的两种方式2.