10-19：数字信号处理课程小结

怕什么真理无穷，进一寸有一寸的欢喜。 ——胡适

1.再论频率分辨率

以下代码均为课上钱老师代码，放着收藏。

import matplotlib.pyplot as plt
import numpy as np

#ground truth spec
#用100s的信号来近似代表真实的谱图
'''
sigL1/fs等于N*时间间隔，算出来就是信号物理时间长度
'''
sigL1 = 100000
fs = 1000
nfft1 = sigL1

f1 = 100
f2 = 100.1
f3 = 20.4

t1 = np.arange(0,sigL1)/fs                     # 采样的时刻
#the point number nfft//2 reprents fs/2
half_fs1 = nfft1//2                            # 这里是取一半，因为FFT是镜像对称的
ff1 = fs*np.arange(0,half_fs1+1)/(2*half_fs1)  # 这里是算频率

sig1 = 0.5*np.sin(2*np.pi*f1*t1)+0.8*np.sin(2*np.pi*f2*t1)+0.7*np.sin(2*np.pi*f3*t1)+0.1*np.random.randn(sigL1)
s1 = np.fft.fft(sig1,nfft1)/sigL1              # 做FFT变换，并归一化
p1 = 10*np.log10(np.real(s1)**2+np.imag(s1)**2)# 算dB，如果log里面没有平方，前面系数就是20，有平方就是10
plt.figure(figsize=(12,4))
plt.plot(ff1,p1[:half_fs1+1],linewidth=2,label="simulated ground truth")
plt.legend()

#看一下20Hz附近
seg1 = np.arange(2000,2100)
plt.figure(figsize=(12,4))
plt.plot(ff1[seg1],p1[seg1],linewidth=2,label="simulated ground truth")
plt.legend()
#看一下100Hz附近
seg2 = np.arange(9950,10050)
plt.figure(figsize=(12,4))
plt.plot(ff1[seg2],p1[seg2],linewidth=2,label="simulated ground truth")
plt.legend()

################################################################################
#信号长度1000点，fs=1000，fft也是1000点，即1s完整信号
sigL2 = 1000
nfft2 = sigL2
t2 = np.arange(0,sigL2)/fs                      # 同样的，采样时刻
half_fs2 = nfft2//2
ff2 = fs*np.arange(0,half_fs2+1)/(2*half_fs2)   # 同样计算频率
sig2 = 0.5*np.sin(2*np.pi*f1*t2)+0.8*np.sin(2*np.pi*f2*t2)+0.7*np.sin(2*np.pi*f3*t2)+0.1*np.random.randn(sigL2)
s2 = np.fft.fft(sig2,nfft2)/sigL2
p2 = 10*np.log10(np.real(s2)**2+np.imag(s2)**2)

plt.figure(figsize=(12,4))
plt.plot(ff2,p2[:half_fs2+1],linewidth=2,label="1s,fs=1kHz")
plt.legend()

#看一下100Hz附近，不仅如此，还将上一个fft的结果放在一起比较
plt.figure(figsize=(12,4))
plt.plot(ff1[9800:10200],p1[9800:10200],linewidth=2,label="simulated ground truth")
plt.plot(ff2[98:103],p2[98:103],'ro',label="1s,fs=1kHz")
plt.legend()
#看一下20Hz附近
plt.figure(figsize=(12,4))
plt.plot(ff1[1800:2100],p1[1800:2100],linewidth=2,label="simulated ground truth")
plt.plot(ff2[18:22],p2[18:22],'ro',label="1s,fs=1kHz")
plt.legend()

################################################################################
#信号长度1000点，fs=1000，fft是100000点，即补零到100s
# sig2是1000点的信号，补到100000，但是归一化除于1000，如果只是归一化的话，除于100000应该也可以
s3 = np.fft.fft(sig2,nfft1)/sigL2               
p3 = 10*np.log10(np.real(s3)**2+np.imag(s3)**2)

plt.figure(figsize=(12,4))
plt.plot(ff1,p3[:half_fs1+1],linewidth=2,label="1s,fs=1kHz,1e5 nfft")
plt.legend()
#看一下20Hz附近
plt.figure(figsize=(12,4))
plt.plot(ff1[seg1],p1[seg1],linewidth=2,label="simulated ground truth")
plt.plot(ff1[seg1],p3[seg1],'ro',label="1s,fs=1kHz,1e5 nfft")
plt.legend()
#看一下100Hz附近
plt.figure(figsize=(12,4))
plt.plot(ff1[seg2],p1[seg2],linewidth=2,label="simulated ground truth")
plt.plot(ff1[seg2],p3[seg2],'ro',label="1s,fs=1kHz,1e5 nfft")
plt.legend()

老师取频率为100,100.1和20.4的正弦信号和随机噪声叠加成的信号进行分析，看看取不同的物理时间长度和补0能不能分辨出这三个频率。结论是：1.补0并不能影响物理分辨率，调采样频率也没用，物理分辨率只跟信号物理时间长度有关，比如100s VS 1s。2.补0无法提升物理分辨率，但是提升了横轴的计算精度！其实是相当于在结果间隙进行插值。

2.时域处理

这里同样是老师上课讲的代码，我增加了一些对自己有帮助的注释。

import numpy as np
import matplotlib.pyplot as plt
from numpy.lib.stride_tricks import as_strided

#读取实验数据并且转换成一维数组
sig = np.load("LabIsig.npy")
sig = sig.squeeze()        # shape中为1的维度去掉

#数据N点
N = 4000
#采样率2000Hz
fs = 2000
#时间标签
t = np.arange(0,N)*(1/fs)  # 采样的时刻


#绘制信号
plt.figure(figsize=(10,5))
plt.plot(t,sig,linewidth=2,label="signal")

#坐标标注和调整
plt.xlabel("Time(s)")
plt.ylabel("Volt")
plt.ylim(-0.5,0.5)
#显示图例
plt.legend()
plt.show()

#归一化到1
sig = sig/np.max(sig)

#信号分段
#计算需要多少段
segL = 64               # 这个是每一小段的长度
overlap = 10            
delta = segL-overlap

# 这里算需要多少段，(N-overlap)/(M-overlap ),M表示段长
segNum = np.int32(np.ceil((N-overlap)/delta));

#扩展信号：看最后有没有多出来一点，补0处理
padNum = segNum*delta+overlap-N
if padNum==0:
    sigEx = sig
elif padNum>0:
    sigEx = np.hstack((sig,np.zeros(padNum)))    

#分段标签:其实就是找每一段的起始位置
segIdx = np.arange(0,segNum)*delta
#生成分段矩阵
'''
1.关于array.strides的解释：https://numpy.org/doc/stable/reference/generated/numpy.ndarray.strides.html
  貌似在这里是每过8个字节可以移动到下一个位置
2.关于as_strided()函数的解释：https://numpy.org/doc/stable/reference/generated/numpy.lib.stride_tricks.as_strided.html#numpy.lib.stride_tricks.as_strided
  在这里sigEx.strides[0]*delta = 432；sigEx.strides[0] = 8；sigEx.shape=(4006,)
  这是在说，我需要移动8字节到下一个数，那我移动8*54个字节，这样就跟上面segIdx对上了。最后得到的segMat大小为segNum*segL（74*64）
  由此完成分段，得到一个矩阵
'''
segMat = as_strided(sigEx,shape=(segNum,segL),strides=(sigEx.strides[0]*delta,sigEx.strides[0]))

# 局部平均能量
E1 = np.sum(segMat**2,axis=1)/segL
# 局部最大值
E2 = np.max(segMat,axis=1)
'''
axis = 0,我们总是习惯先说行再说列，axis=0即我们把一行看成一个运算单元，得到的结果也是相同大小的一行
axis = 1，我们把一列看成运算单元，得到的结果跟一列大小相同
'''
plt.figure()
plt.plot(E1)
plt.plot(E2)
plt.legend(['E1','E2'])
plt.show()

loc_pow = E1
loc_max = E2
plt.figure(figsize=(10,5))
plt.plot(t,sig, linewidth=2, label="signal without noise")
# t[segIdx]取法就是74段开始时所在的时刻，目的是为了将图缩放在一起
plt.plot(t[segIdx],loc_pow, color="red", linewidth=2, label="local mean power")
plt.plot(t[segIdx],loc_max,"r--",linewidth=2,label="local max")
plt.legend()
plt.show()

pow_diff=np.convolve(loc_pow,np.array([1,-1]),'full')
pIdx=[]
for i in range(segNum):
    if pow_diff[i]* pow_diff[i+1]<0 and loc_pow[i]>0.05:
        pIdx.append(i)
plt.plot(loc_pow)
plt.scatter(pIdx,loc_pow[pIdx],c='r')
plt.show()

这里处理的信号是一个心音信号，如下图展示。想要区分S1和S2。这一小节叫时域处理，看来没频率啥啥事了，让我们看看只处理时域能不能区分S1和S2。

心音信号.png

老师教我们这样处理：先进行分段，然后计算每一段的局部平均能量

E_1=\sum x^2

和每一段的局部最大值

E_2=max |x|

，

x

指一段中的时域信号。然后我们找

E_1

的峰值，S1和S2特征会出现峰值的差异。由此我们大致区分开了S1和S2。

时域处理后的结果.png

3.STFT

以下代码大部分来自课上老师：

# -*- coding: utf-8 -*-
"""
Created on Sun Oct 18 16:22:25 2020

@author: qx-HW
"""
import matplotlib.pyplot as plt
import numpy as np
from numpy.lib.stride_tricks import as_strided
#时频原子参数
N = 3000
fs = 1000
alpha = 100
t = np.arange(0,N)*(1/fs)
t0 = 0.4
f0 = 0
f1 = 100
f = (f1-f0)*t+f0

sig = ((alpha/np.pi)**0.25)*np.exp(-0.5*alpha*(t-t0)*(t-t0))*np.sin(2*np.pi*f*t)
plt.figure(figsize=(12,4))
plt.plot(t,sig,linewidth=2)
plt.show()

#fft点数以及频率轴
nfft = N
half_fs = nfft//2
ff = fs*np.arange(0,half_fs+1)/(2*half_fs)

segL = 32               # 这个是每一小段的长度
overlap = 10     
def FenDuan(segL,overlap,sig,N):
    #信号分段
    #计算需要多少段 
    delta = segL-overlap
    
    # 这里算需要多少段，(N-overlap)/(M-overlap ),M表示段长
    segNum = np.int32(np.ceil((N-overlap)/delta));
    
    #扩展信号：看最后有没有多出来一点，补0处理
    padNum = segNum*delta+overlap-N
    if padNum==0:
        sigEx = sig
    elif padNum>0:
        sigEx = np.hstack((sig,np.zeros(padNum)))    
    
    #分段标签:其实就是找每一段的起始位置
    segIdx = np.arange(0,segNum)*delta
    #生成分段矩阵
    segMat = as_strided(sigEx,shape=(segNum,segL),strides=(sigEx.strides[0]*delta,sigEx.strides[0]))
    return segMat,segIdx

segMat,segIdx = FenDuan(segL,overlap,sig,N)

pMat = np.fft.fft(segMat) # (19,64)

'''
coutour([X, Y,] Z,[levels], **kwargs)
https://matplotlib.org/api/_as_gen/matplotlib.pyplot.contourf.html
是来绘制等高线的函数，X,Y确定低下的平面坐标，Z确定三维的高度坐标，levels确定轮廓线/区域的数量和位置
'''
plt.contourf(t[segIdx], ff[:pMat.shape[1]], pMat[:,0:half_fs+1].T,50)
plt.show()

plt.contourf(t[segIdx], ff[:pMat.shape[1]//2], pMat[:,:pMat.shape[1]//2].T,30)
plt.show()

image.png

从代码上来看，就是对一个信号进行分段，然后每段做FFT，然后可视化出来。

集合网络资料和课本对短时傅里叶变换进行一些梳理和总结：

通过时间窗内的一段信号来表示某一时刻的信号特征。窗的长度决定频谱图的时间分辨率和频率分辨率。窗长越长，截取的信号越长。
信号越长，傅里叶变换后的频率分辨率越高，时间分辨率越差。相反，窗长越短，截取的信号就越短，频率分辨率越差，时间分辨率越好。
时间和频率是一对不可兼得的矛盾体，所以绝对意义的瞬时频率是不存在的，只能分时段。在短时傅里叶变换中，时间分辨率和频率分辨率之间不能兼得，应该根据具体需求进行取舍。

短时傅里叶变换就是先把一个函数和窗函数进行相乘，然后再进行一维的傅里叶变换。并通过窗函数的滑动得到一系列的傅里叶变换结果，将这些结果竖着排开得到一个二维的表象。

傅里叶变换后，横轴为频率，纵轴为幅值。短时傅里叶变换后，横轴为时间，纵轴为频率