按时间抽选的fft算法golang实现

按时间抽选的fft算法golang实现

参考博客：https://zhuanlan.zhihu.com/p/135259438

实现目的：参考以上链接博客c语言的按时间抽选的fft算法进行学习，用golang进行重构简化代码并且完成实现功能。

代码：

package main

import (
	"fmt"
	"math"
	"os"
	"strconv"
)

const (
	PI float64 = 3.1415926535
	N          = 128
	M          = 7
	A0 float64 = 255.0
)

var SinIn [N]float64 //正弦波
var DataR [N]float64 //输入波形实数部分
var DataI [N]float64 //输入波形虚数部分
var DataA [N]float64 //输出频谱

func main() {
	InputWave()   //写入输入波形,打点数据存储在InputWave.txt
	ReverseBits() //码元倒序0-127
	FFT()         //FFT计算
	OutputFre()   //输出频谱打点输出至OutputFre.txt
}

//写入输入波形
func InputWave() {
	for i := 0; i < N; i++ {
		SinIn[i] = A0 * (math.Sin(2*PI*float64(i)/25.0) + math.Sin(2*PI*float64(i)*0.4))
		DataR[i] = SinIn[i]
		DataI[i] = 0
		DataA[i] = 0

	}
	//输入波形打点写入InputWave.txt
	file, err := os.OpenFile("InputWave.txt", os.O_CREATE|os.O_TRUNC|os.O_WRONLY, 0666)
	if err != nil {
		fmt.Println("open file failed, err:", err)
		return
	}
	defer file.Close()
	for i := 0; i < N; i++ {
		DataR := strconv.FormatFloat(DataR[i], 'f', 8, 64) //将float64--->string
		file.Write([]byte(DataR + "\n"))
	}
}

//码元倒序0-127
func ReverseBits() {
	var I, J, F0, F1, m, n uint
	for I = 0; I < N; I++ { //对数组元素执行码间倒序

		//获取下标I的反序J的数值
		J = 0
		var k int

		for k = 0; k < (M+1)/2; k++ { //k表示操作
			//反序操作
			m = 1                      //m是最低位为1的二进制数
			n = uint(math.Pow(2, M-1)) //n是第M位为1的二进制数
			m <<= k                    //对m左移k位
			n >>= k                    //对n右移k位
			F0 = I & n                 //I与n按位与提取出前半部分第k位
			F1 = I & m                 //I与m按位与提取出F0对应的后半部分的低位
			if F0 != 0 {
				J = J | m //J与m按位或使F0对应低位为1
			}
			if F1 != 0 {
				J = J | n //J与n按位或使F1对应高位为1
			}
		}
		fmt.Printf("I为：")
		PrintBin(I, M)
		fmt.Printf(" %d\n", I)
		fmt.Printf("J为：")
		PrintBin(J, M)
		fmt.Printf(" %d\n\n", J)
		var temp float64
		if I < J {
			temp = DataR[I]
			DataR[I] = DataR[J]
			DataR[J] = temp
		}
	}
}

//打印二进制串 i为10进制数,bits为二进制位数
func PrintBin(i uint, bits int) {
	for bits != 0 {
		bits--
		if i&(1<<bits) != 0 {
			fmt.Print("1")
		} else {
			fmt.Print("0")
		}

	}
}

//FFT运算

func FFT() {
	var Tr, Ti float64

	for L := 1; L <= M; L++ {
		//1-M次蝶形运算
		//第L级蝶形运算,每个分组个数B为2^(L-1)
		//旋转指数p
		//旋转指数p增量k=2^(M-L)
		var k int
		var B int
		var p float64
		var r int
		B = int(math.Pow(2, float64(L-1)))
		for j := 0; j <= B-1; j++ {
			k = int(math.Pow(2, float64(M-L)))
			p = float64(j * k)
			for i := 0; i <= k-1; i++ {
				//数组下标定为r
				r = j + 2*B*i
				Tr = DataR[r+B]*math.Cos(2.0*PI*p/float64(N)) + DataI[r+B]*math.Sin(2.0*PI*p/float64(N))
				Ti = DataI[r+B]*math.Cos(2.0*PI*p/float64(N)) - DataR[r+B]*math.Sin(2.0*PI*p/float64(N))
				DataR[r+B] = DataR[r] - Tr
				DataI[r+B] = DataI[r] - Ti
				DataR[r] = DataR[r] + Tr
				DataI[r] = DataI[r] + Ti
			}
		}
	}

}

//计算存入幅度谱
func OutputFre() {

	for i := 0; i < N; i++ {
		DataA[i] = math.Sqrt(DataR[i]*DataR[i] + DataI[i]*DataI[i])
	}
	//输出频谱打点至 OutputFre.txt
	file, err := os.OpenFile("OutputFre.txt", os.O_CREATE|os.O_TRUNC|os.O_WRONLY, 0666)
	if err != nil {
		fmt.Println("open file failed, err:", err)
		return
	}
	defer file.Close()
	for i := 0; i < N; i++ {
		DataA := strconv.FormatFloat(DataA[i], 'f', 8, 64) //将float64--->string
		file.Write([]byte(DataA + "\n"))
	}

}

其中最为复杂是蝶形运算中的状态转移方程（迭代方程）：根据书上

但是旋转因子不好直接表示，需要转化一下方程：

（旋转因子P表示为每一级的增量k * j ，j是遍历每个分组个数B，0-B-1）

		Tr = Xr(J+B)cos(2.0*PI*p/N) + Xi(J+B)sin(2.0*PI*p/N)
	    Ti = Xi(J+B)cos(2.0*PI*p/N) - Xr(J+B)sin(2.0*PI*p/N)
		Ar[J] = Xr(J) + Tr 
	    Ai[J] = Xi(J) + Ti
	    Ar[J+B] = Xr(J) - Tr
	    Ai[J+B] = Xi(J) - Ti
	   (其中 Xr为上一级的Ar, Xi为上一级的Ai)

最复杂的数学推导复制一下博客内容：

利用打点数据输出波形在matlab中显示：