用python实现SDFT

上一篇文章介绍了如何用python实现DFT，本篇文章采用python实现SDFT。
首先，生成一个由180Hz，390Hz和600Hz的正弦波叠加的采样信号，幅值分别为7.0，2.8，5.1，采样时间设置为1s，采样频率设置为4000Hz。

Ts = 1
fs = 4000
xs = np.linspace(0, Ts, int(Ts*fs))  # 生成时间
# 生成采样信号 由180Hz，390Hz和600Hz的正弦波叠加
ys = 7.0*np.sin(2*np.pi*180*xs) + 2.8*np.sin(2*np.pi*390*xs) + 5.1*np.sin(2*np.pi*600*xs)
plt.plot(xs, ys)
plt.show()

信号时域图如下

定义自己的SDFT类：

class mySDFT:
    def __init__(self, k, M, ys, fs):
        '''
            :param k: 频域索引
            :param M: 序列长度
            :param ys: 采样信号
            :param fs: 采样频率
            :return:
        '''
        self.k = k
        self.M = M
        self.ys = ys
        self.fs = fs

    def cal_Wm(self, p):
        return np.cos(2*np.pi*p/self.M) + 1j*np.sin(2*np.pi*p/self.M)

    def DFT(self, date):
        Xnk_pre = []
        for k in range(self.k):
            # 计算Xnk
            Xnk = 0
            for m in range(self.M):
                Xnk += date[m] * self.cal_Wm(-(m*k))
            Xnk_pre.append(Xnk)  # shape:[self.k, 1]

        return Xnk_pre

    def SDFT(self):
        date = self.ys[:self.M]
        Xnk_pre = self.DFT(date)

        res = []
        Xnk_new = []

        for i in range(1, len(self.ys)-self.M):
            date = self.ys[i:self.M + i]
            for k in range(self.k):
                Xnk_new.append(self.cal_Wm(k) * (Xnk_pre[k] - date[0] + date[-1]))

            Xnk_pre = copy.copy(Xnk_new)
            Xnk_new.clear()

        res.append(Xnk_pre)
        res = np.array(res)
        res = res.reshape(-1, self.k)

        return res

    def plotPPT(self, data):
        # 绘制频谱图
        A = abs(np.array(data))  # 计算幅度

        amp_x = A / self.M * 2  # 纵坐标变换
        label_x = np.linspace(0, int(self.M / 2) - 1, int(self.M / 2))  # 生成频率坐标
        amp = amp_x[0:int(self.M / 2)]  # 选取前半段计算结果即可，幅值  对称

        fs = self.fs  # 计算采样频率
        fre = label_x / self.M * fs  # 频率坐标变换
        
        # 存储结果
        save_data = {'frequence':fre, 'amp':amp}
        df = pd.DataFrame(save_data)
        path = './Data/data_SDFT.csv'
        df.to_csv(path, index=False)
        
        # 绘制频谱图
        plt.figure()
        plt.plot(fre, amp)
        plt.title('Amplitute-Frequence-Curve')
        plt.ylabel('Amplitute / a.u.')
        plt.xlabel('Frequence / Hz')
        plt.show()

其中的SDFT的计算公式为：

实例化自己定义的SDFT类，并将初始的M值设置为400，即取前400个点进行DFT计算，后面的3600点由上面的公式迭代算出。

M = 400
k = M
sdft = mySDFT(k,M,ys,fs)
res = sdft.SDFT()
sdft.plotPPT(res[-1, :])

绘制结果频谱图：

从结果可以看出，频率大致可以分辨出来，但幅值差别较大。
尝试将M设置为2000，绘制结果频谱图：

此时结果与真实值已非常接近，将此次结果进行傅里叶反变换：

plt.plot(np.arange(len(res[-1, :])), ifft(res[-1, :]))
plt.show()

初探傅里叶变换，如有不足之处，望不吝赐教。