dft code 

 1 #!/usr/bin/env python
 2 # -*- encoding: utf-8 -*-
 3 '''
 4 @File    :   dft_test.py
 5 @Time    :   2023/08/02 08:49:41
 6 @Author  :   qfr
 7 @Version :   1.0
 8 @Contact :   qfr@whu.edu.cn
 9 @License :   (C)Copyright 2021-2025, qfr&lz
10 @Desc    :   None
11 '''
12 
13 # here put the import lib
14 import numpy as np
15 import matplotlib.pyplot as plt
16 
17 def GenSignal(CN0, sigMa, coh, N):
18     snr = np.power(10, CN0 / 10) / (1 / coh)
19     a = np.sqrt(snr * 2 * sigMa ** 2)
20 
21     data = np.zeros((N,), dtype=np.complex128)
22     data.real = sigMa * np.random.randn(N)
23     data.imag = sigMa * np.random.randn(N)
24     data += a * np.ones((N,))
25     return data
26 
27 def DftTest(CN0, sigMa, coh, N):
28     fs = 1 / coh
29     data = GenSignal(CN0, sigMa, coh, N)
30 
31     # test type 1
32     dftFreqNum = N                            # try to modify N
33     dftPhaseResolution = 2 * np.pi / N        # try to modify N
34     '''
35     # test type 2
36     dftFreqNum = 1000                            # try to modify N
37     dftPhaseResolution = 2 * np.pi / 1000        # try to modify N
38     '''
39     dftData = []
40     dftFreq = []
41     for i in range(dftFreqNum):
42         deltPhae = dftPhaseResolution * i
43         # accord phase to calculate dft Freq
44         if deltPhae < np.pi:  # positive freq
45             tmpFreq = deltPhae / (2 * np.pi) * fs
46             dftFreq.append(tmpFreq)
47         elif deltPhae < 2 * np.pi:  # negative freq
48             tmpFreq = (deltPhae - 2 * np.pi) / (2 * np.pi) * fs
49             dftFreq.append(tmpFreq)
50         else:
51             print("current para is error, not fit [dftFreqNum * dftPhaseResolution <= 1]")
52             continue
53 
54         tmpData = 0
55         tmpPhase = 0
56         for d in data:
57             tmpPhase += deltPhae
58             tmp = np.cos(tmpPhase) - 1j * np.sin(tmpPhase)
59             tmpData += d * tmp
60 
61         dftData.append(tmpData)
62 
63     #  dft频率大小排列，为绘图展示
64     #  原因：和fft一样，变换的结果是正频率在前，负频率在后，调整恢复从大到小的顺序
65     dftFreq = np.array(dftFreq)
66     dftData = np.array(dftData)
67     sortIdx = np.argsort(dftFreq)
68     dftFreq = dftFreq[sortIdx]
69     dftData = dftData[sortIdx]
70 
71     # 绘图
72     fig = plt.figure(1)
73     ax = fig.subplots(nrows=2, ncols=1, sharex='col')
74     ax1 = ax[0]     # abs
75     ax2 = ax[1]     # phase
76 
77     ax1.plot(dftFreq, np.abs(dftData),c='red', label='abs',linewidth = 1,linestyle='-', marker = '*')#c:颜色 label:图例
78     ax2.plot(dftFreq, np.rad2deg(np.angle(dftData)), c='red', label='phase',linewidth = 1,linestyle='-', marker = '*')#c:颜色 label:图例
79 
80     ax1.set_title("dft result cn0:" + str(CN0),size=22)
81     ax1.grid()
82     ax1.set_ylabel("mag ",size=22)
83     ax1.legend(loc='upper left')
84 
85     ax2.legend(loc='upper left')
86     ax2.grid()
87     ax2.set_ylabel("phase /deg",size=22)
88     ax2.set_xlabel("freq /Hz",size=22)
89     plt.show()
90 
91 
92 if __name__ == '__main__':
93     DftTest(35, 1, 0.02, 5)

 

posted @ 2023-08-02 09:53  博客园—哆啦A梦  阅读(36)  评论(0)    收藏  举报