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

参照destiny0321的《合成孔径雷达成像——算法与实现》之【1】仿真图2.2给出的Matlab程序，用Python试了一下。

import numpy as np
import math
import cmath
import pylab
pylab.mpl.rcParams['font.sans-serif'] = ['SimHei']
pylab.mpl.rcParams['axes.unicode_minus']=False    #用来正常显示负号
## 参数设置
M = 256                 # 矩阵高度
N = 256                 # 矩阵宽度
top = M/8+1 
bottom = M*7/8 
left = N/8+1 
right = N*7/8 
theta = math.pi/12           # 扭曲或旋转角度

## 生成信号
# 原始信号
S0 = np.zeros((M,N))
S0[int(top):int(bottom),int(left):int(right)] = 1 
#S0[256/8+1:256*7/8 ,256/8+1:256*7/8] = 1 
# 扭曲信号
S1 = np.zeros((M,N)) 
for ii in range(0,M):
    for jj in range(0,N):
        x = jj-N/2 
        y = (M+1-ii)-M/2 
        xx = round(x+N/2)
        yy = M+1-round(x*math.sin(-theta)+y*math.cos(-theta)+M/2)
        if(yy>=0 and yy<= M-1):
            S1[ii,jj] = S0[yy,xx]
# 旋转信号
S2 = np.zeros((M,N)) 
for ii in range(0,M):
    for jj in range(0,N):
        x = jj-N/2 
        y = (M+1-ii)-M/2 
        xx = round(x*math.cos(-theta)-y*math.sin(-theta)+N/2) 
        yy = M+1-round(x*math.sin(-theta)+y*math.cos(-theta)+M/2) 
        if(xx>=0 and xx<= N-1 and yy>=0 and yy<=M-1):
            S2[ii,jj] = S0[yy,xx] 
## 二维傅里叶变换
# 原始信号的二维傅里叶变换
S0_ff = np.fft.fftshift(np.fft.fft2(np.fft.fftshift(S0))) 
S0_ff = abs(S0_ff) 
S0_ff = S0_ff/S0_ff.max() 
S0_ff = 20*np.log10(S0_ff+1e-4) 
# 原始信号二维傅里叶变换
S1_ff = np.fft.fftshift(np.fft.fft2(np.fft.fftshift(S1))) 
S1_ff = abs(S1_ff) 
S1_ff = S1_ff/S1_ff.max() 
S1_ff = 20*np.log10(S1_ff+1e-4) 
# 原始信号二维傅里叶变换
S2_ff = np.fft.fftshift(np.fft.fft2(np.fft.fftshift(S2))) 
S2_ff = abs(S2_ff) 
S2_ff = S2_ff/S2_ff.max() 
S2_ff = 20*np.log10(S2_ff+1e-4) 

## 画图
pylab.figure(1)
pylab.plt.subplot(2,3,1)
pylab.pcolormesh(S0,cmap='jet')
pylab.plt.gca().invert_yaxis()
pylab.title('（a）时域，原始信号')
pylab.plt.subplot(2,3,4)
pylab.pcolormesh(S0_ff,cmap='jet')
pylab.plt.gca().invert_yaxis()
pylab.title('（b）原始信号频谱')
pylab.plt.subplot(2,3,2)
pylab.pcolormesh(S1,cmap='jet')
pylab.plt.gca().invert_yaxis()
pylab.title('（c）时域，扭曲信号')
pylab.plt.subplot(2,3,5)
pylab.pcolormesh(S1_ff,cmap='jet')
pylab.plt.gca().invert_yaxis()
pylab.title('（d）扭曲信号频谱')
pylab.plt.subplot(2,3,3)
pylab.pcolormesh(S2,cmap='jet')
pylab.plt.gca().invert_yaxis()
pylab.title('（e）时域，旋转信号')
pylab.plt.subplot(2,3,6)
pylab.pcolormesh(S2_ff,cmap='jet')
pylab.plt.gca().invert_yaxis()
pylab.title('（f）旋转信号频谱')
pylab.show()