FRI信号重构算法（一）：零化滤波器法（python版）

## FRI 信号重构
## 信号模型为Dirac脉冲，采样核为sinc核，直接取傅里叶系数
## 
import numpy as np
import matplotlib.pyplot as plt
import math
fs = 1000 #采样率
dt = 1/fs #采样间隔
simuT = 1 #仿真时间
L = 5 #Dirac脉冲数，2*L即新息率
simuTime = np.arange(0,simuT,dt) #仿真时间向量，步长为dt，起点为0，终点为1，不包含终点
ak = np.array([0.4, 0.6, 0.3, 0.9, 0.7])  #幅值
tk = np.array([0.1, 0.2, 0.3, 0.7, 0.8]) #时延,单位为秒
tk_index = (tk*fs).astype(np.int32) #脉冲流信号的脉冲位置的下标，下标必须是整数，astype()用于强制类型转换
N = simuTime.size  #信号的长度
x = np.zeros((N))  #生成全零的仿真信号
x[tk_index] = ak   #添加脉冲
## 绘制仿真信号
#plt.stem(simuTime,x)
#plt.show()
X = np.fft.fftshift(np.fft.fft(x))/N  #fft变换，并且将零分量平移到中心

mid = (np.ceil(N/2)).astype(int) #零频率分量的下标
X_part = X[mid:mid+2*L] #采用sinc核，频域为理想低通滤波器，因此只取中间的部分的傅里叶系数即可
# toeplitz矩阵生成函数
# 输入参数: c : 第一列
#          r :  第一行
# 返回参数  T ：toeplitz矩阵
def toeplitz(c,r):
    m = c.size #toeplitz矩阵的行数
    n = r.size #toeplitz矩阵的列数
    if c.dtype == complex:
        T = np.ones([m,n],dtype=complex)
    else:
        T = np.ones([m,n])
    c_rev = c[::-1] #翻转
    c_rev_r = np.zeros(m+n-1,dtype=complex) #拼接向量,因为首元素相同，所以拼接后向量的长度为m+n-1
    c_rev_r[0:m] = c_rev
    c_rev_r[m:m+n-1] = r[1:]
    for i in range(m):
        T[i,:] = c_rev_r[m-i-1:m-i-1+n]
    return T
c = X_part[(L-1):(2*L-1)] #toeplitz矩阵的第一列
r = X_part[(L-1)::-1]     #toeplitz矩阵的第一行
T = toeplitz(c,r) #获得toeplitz矩阵
## 求解yelu—walker方程 TA=y
y = X_part[L:]
y = y.reshape([L,1])
A = np.dot((np.linalg.inv(T)),-y) #注意矩阵求逆及矩阵乘法
A1 = np.zeros([A.size+1,1],dtype=A.dtype) #给向量添加一个1
A1[0,0] = 1+0j   #添加1
A1[1:,0] = A[:,0] #补充后面的系数
z = np.roots(A1[:,0]) #求根
tkr = np.sort(np.mod((-np.angle(z)*simuT/(2*math.pi)),simuT))

print(tkr)
plt.stem(tk,ak,markerfmt='D')
plt.stem(tkr,ak,markerfmt='*')

plt.show()