Python实现《合成孔径雷达成像——算法与实现》图3.8

Python实现《合成孔径雷达成像——算法与实现》图3.8。

import matplotlib.pyplot as plt
import numpy as np
import math
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus']=False    #用来正常显示负号
#信号持续时间T=7.24us，信号带宽B=5.8MHz，将过采样率设为5是为了更清晰地观测信号波形
T = 7.24e-6                 # 信号持续时间
B = 5.8e6                   # 信号带宽
K = B/T                     # 调频率
ratio = 10                  # 过采样率
Fs = ratio*B                # 采样频率
dt = 1/Fs                   # 采样间隔
N = math.ceil(T/dt)              # 采样点数
tc = 1e-6                   # 脉冲中心相当于t=0的时间偏移
t = np.arange((0-N/2)/N*T,(N-N/2)/N*T,dt) # 时间轴
st = np.exp(1j*math.pi*K*np.multiply(t-tc, t-tc))      #生成信号
ht = np.exp(-1j*math.pi*K*np.multiply(-t-tc, -t-tc))      # 匹配滤波器为什么-tc不是+tc?
out = np.fft.fftshift(np.fft.ifft(np.multiply(np.fft.fft(st),np.fft.fft(ht))))
Z = np.abs(out)
Z = Z/max(Z)
Z = 20*np.log10(np.finfo(np.float64).eps+Z)


## 画图
plt.figure(1)
plt.subplot(2,2,1)
plt.title('（a）输入阵列信号的实部')
plt.plot(t*1e6,np.real(st))
plt.ylabel('幅度')

plt.subplot(2,2,2)
plt.title('(c)压缩后的信号(经扩展)')
plt.axis([-1,1,-30,0])
plt.plot(t*1e6,Z)
plt.ylabel('幅度(dB)')

plt.subplot(2,2,3)
plt.title('（b）压缩后的信号')
plt.plot(t*1e6,out)
plt.xlabel('相对于$\mathregular{t_0}$''时间(μs)')
plt.ylabel('幅度')

plt.subplot(2,2,4)
plt.title('（d）压缩后信号的相位(经扩展)')
plt.axis([-1,1,-5,5])
plt.plot(t*1e6,np.angle(out))
plt.xlabel('相对于$\mathregular{t_0}$''时间(μs)')
plt.ylabel('相位（弧度）')

plt.show()

#加窗减小旁瓣
window=np.kaiser(N, 2.5)#典型值β=2.5见书P60
ht_window = np.multiply(window,ht)   #时域加窗匹配滤波器

Sf =  np.fft.fftshift( np.fft.fft( np.fft.fftshift(st)))
Hf =  np.fft.fftshift( np.fft.fft( np.fft.fftshift(ht)))

Hf_window = np.fft.fftshift(np.fft.fft(np.fft.fftshift(ht_window)))

out = np.fft.ifftshift(np.fft.ifft(np.fft.ifftshift(np.multiply(Sf,Hf))))
out_window = np.fft.ifftshift(np.fft.ifft(np.fft.ifftshift(np.multiply(Sf,Hf_window))))


Z1 = np.abs(out)
Z1 = Z1/max(Z1)
Z1 = 20*np.log10(Z1)

Z2 = np.abs(out_window)
Z2 = Z2/max(Z2)
Z2 = 20*np.log10(Z2)

tt = np.linspace(-0.5,0.5,N)
plt.figure(2)
plt.subplot(2,2,1)
plt.title('脉冲压缩之后的信号（未加窗）')
plt.plot(tt,out)
plt.ylabel('幅度')

plt.subplot(2,2,2)
plt.title('脉冲压缩之后的信号（未加窗）')
plt.axis([-0.3,0.3,-35,0])
plt.plot(tt,Z1)
plt.ylabel('幅度(dB)')

plt.subplot(2,2,3)
plt.title('脉冲压缩之后的信号（加窗）')
plt.plot(tt,out_window)
plt.xlabel('时间（归一化后）')
plt.ylabel('幅度')

plt.subplot(2,2,4)
plt.title('脉冲压缩之后的信号（加窗）')
plt.axis([-0.3,0.3,-35,5])
plt.plot(tt,Z2)
plt.xlabel('时间（归一化后）')
plt.ylabel('幅度(dB)')

plt.show()