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

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

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 = 42e-6                 # 信号持续时间见书P61
B = 17.2e6                  # 信号带宽
K = B/T                     # 调频率
ratio = 1.07                # 过采样率
Fs = ratio*B                # 采样频率
dt = 1/Fs                   # 采样间隔
N = math.ceil(T/dt)         # 采样点数
print(N)
t = np.arange((0-N/2)/N*T,(N-N/2)/N*T,dt)       # 时间轴
st0 = np.exp(1j*math.pi*K*np.multiply(t, t))      # 生成基本信号

N1 =len(t)
print(N1)

f = np.linspace((0-N1/2)/N1*Fs,(N1-N1/2)/N1*Fs,N1)     # 频率轴
print(len(f))
Sf = np.fft.fftshift(np.fft.fft(st0))                # FFT

Nfft=2048
window = np.kaiser(N1,2.5)    # 时域窗
space = np.zeros(Nfft-N1)
a = window*st0
b = space
e = [a,b]#补0
e = [y for x in e for y in x]
#e = window*e
Hf2 = np.conj(np.fft.fft(e))   # 方式2的匹配滤波器：补零后计算DFT，对结果取复共轭
#Hf2 = np.fft.fftshift(np.conj(np.fft.fft(e)))

window = np.kaiser(N1,2.5)    # 时域窗
Window = np.fft.fftshift(window)  # 频域窗
Hf3 = Window*np.exp(1j*math.pi*np.multiply(f,f)/K)     # 方式3频域匹配滤波器

Out3 = np.fft.ifft(np.fft.ifftshift(Sf*Hf3))

instantaneous_frequency = np.diff(np.unwrap(np.angle(Hf3)))/dt/(2.0*np.pi) 
plt.figure(1)
plt.subplot(2,1,1)
plt.plot(np.abs(Hf2))
plt.axis('tight')
plt.title('(a)加权后匹配滤波器的幅度谱')
plt.ylabel('幅度')
plt.subplot(2,1,2)
plt.plot(np.unwrap(np.angle(Hf2)))#一般在计算一个系统相频特性时，就要用到反正切函数提取相位，
#计算机中反正切函数规定：在一、二象限中的角度为0 ~ π，三、四象限的角度为0 ~ -π
#但实际得到的结果会发生相位跳变，跳变幅度为2π，
#这就叫做相位的卷绕unwrap的作用就是解卷绕，使相位在π处不发生跳变，从而反映出真实的相位变化
#plt.plot(np.unwrap(np.arctan(np.imag(Hf2),np.real(Hf2) ),2*math.pi))
#plt.plot(np.angle(Hf2))
plt.axis('tight')
plt.title('(b)匹配滤波器的相位谱')
plt.ylabel('相位（弧度）')

plt.figure(2)
plt.subplot(4,1,1)
plt.plot(np.abs(Hf3))
plt.axis('tight')
plt.title('(a)频域匹配滤波器的幅度')
plt.ylabel('幅度')
plt.subplot(4,1,2)
plt.plot(instantaneous_frequency*10e-7)
plt.axis('tight')
plt.title('(b)频域匹配滤波器的瞬时频率')
plt.ylabel('频率（MHz）')
plt.subplot(4,1,3)
plt.plot(np.unwrap(np.angle(Hf3)))
plt.axis('tight')
plt.title('(c)频域匹配滤波器的相位')
plt.ylabel('相位（弧度）')
plt.xlabel('频率(FFT采样点)')
plt.subplot(4,1,4)
plt.plot(np.real(Hf3))
plt.axis('tight')
plt.title('(d)频域匹配滤波器的实部')
plt.xlabel('频率(采样点)')
plt.ylabel('幅度')
plt.show()