圆的FFT分析的算法

ljm1200

已于 2022-09-06 14:13:53 修改

阅读量2k

点赞数

文章标签：算法 python matplotlib pandas

于 2022-09-06 13:56:55 首次发布

本文链接：https://blog.csdn.net/LJM1200/article/details/126723503

版权

傅里叶变换其实是发明了一种基础文字只有两个sin,cos去表达时域信号的语言。也就是说只要懂得sin,cos这两个基础文字，再使用傅里叶变换公式，就能描述所有时域上的信号。

傅里叶级数：

$f(t) = a_0 + \sum _ 0 ^ \infty [a_n cos(nwt) + b_n sin(nwt)]$

对于一个正弦信号可以看成一定频率下的圆周运动，每一时刻下的(X,Y)坐标可以表示为

$X=Rcos(\alpha) Y=Rsin(\alpha)$

所以转换到频域就得到下面的图形

$A_1 sin (2 \pi \times t)$

成分1

频率	幅值	相位
1Hz	A1	0

如果在大圆运动的同时在加一个小圆一起运动：

成分1

频率	幅值	相位
1Hz	A1	0

成分2

频率	幅值	相位
2Hz	A2	0

如果在大圆运动的同时在加一个小圆在小圆上再加一个小圆一起运动：

成分1

频率	幅值	相位
1Hz	A1	0

成分2

频率	幅值	相位
2Hz	A2	0

成分3

频率	幅值	相位
3Hz	A3	0

我们就可以利用这个思想把每段的频率和振幅在频谱图上标示出来。由于采样的信号不是连续的，所以在计算的时候我们采取了“离散时间傅立叶变换”（DTFT）算法，它是一种特定的傅立叶变换，应用于表示接收到的采样信号的离散时间函数x（n）。实际上就是测量圆上点的顺序，和点坐标的离散函数。除了我们为FFT提供N个点样本的这一事实之外，没有太多的差异，DTFT是FT的近似值，必然会在我们的频谱分析中引入精度损失。这被称为“不确定性”

对于一个圆度为0的圆，它的波形成分就应该是一个标准正弦或者余弦波形，由于工件在加工过程中的材料特性，机床精度，刀具，工艺参数等因素的影响，使圆的表面产生了一定的形变，波纹度和粗糙度。那这样就势必会引入其他频率成分的波形，如上图一样有了频率2，频率3，甚至频率n的成分。

对于圆的拟合我们使用最小二乘法（LSQ），这个标准圆是其中的基带成分，那么和这个标准圆的半径上的差值就形成了上面的小圆，所以每个小圆的所引入的频率的叠加，就构成了最终的圆度波形。

在实际中我们测量的是圆上每个点的x，y坐标，这样的函数构不成x(n)的函数，我么通过上面的思想把每个点转换成半径，再计算出半径和标准半径（即最小二乘圆半径）的差值。这样的数据就构成了时间序列和半径差的x(n)的函数，这样我们就可以做傅里叶变换了。

算法如下：假定在360度范围内有M个点，半径为

$R = R_1,R_2,R_3,\cdots R_m$

$\bar R = r_0$

$\Delta r_i = r_i- r_0$

$real_k = \frac{1}{m}(\sum_{i=1}^m(sin(\frac{i\times2\times\pi }{m}\times k)\times \Delta r_i))$

$im_k = \frac{1}{m}(\sum_{i=1}^m(cos(\frac{i\times2\times\pi }{m}\times k)\times \Delta r_i))$

$F_k =\sqrt{ {real_k}^2 + {im_k}^2}$

$(k = 1,2,3, \cdots, n)$

#! /usr/bin/env python
# -*- coding: utf-8 -*-
'''
@Created on: 2022/04/27 19:31:52
@author: ZCJOHNLV
@Version: V1.0
'''
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

# 读取测量点
def readTXT(file):
    f=open(file, encoding='UTF-8')
    df_txt = pd.read_table(f,engine='python',sep=' ',encoding = 'gb18030')
    df = pd.DataFrame(df_txt)
    return df



R0 = 100/2

df = readTXT(r"C:\MySpace\Project\demo\files\123.txt")


x_data = pd.DataFrame(df, columns=['X']).values.tolist()
y_data = pd.DataFrame(df, columns=['Y']).values.tolist()
points_Num = len(x_data)


real = 0
im = 0

AMP = 0

plot_x = []
plot_y = []
    
for j in range(300):
    plot_x.append(j)
    for i in range(points_Num):
        xi = float(x_data[i][0])
        yi = float(y_data[i][0])
        Ri = np.sqrt(xi*xi + yi*yi)
        R_delta = Ri -R0 
        real = real + np.sin((i*2*np.pi)/points_Num*j)*R_delta/points_Num
        im  = im + np.cos((i*2*np.pi)/points_Num*j)*R_delta
    AMP_tmp = np.sqrt(real*real + im*im)
    real = 0
    im = 0
    AMP = 2*round(AMP_tmp,8)
    plot_y.append(AMP)
    print(AMP)

plt.bar(plot_x,plot_y)
plt.show()