在项目中,需要画波形频谱图,因此进行查找,不是很懂相关知识,下列代码主要是针对这篇文章。
http://blog.csdn.net/xcgspring/article/details/4749075
//快速傅里叶变换
/*
入口参数:
inv: =1,傅里叶变换; =-1,逆傅里叶变换
N:输入的点数,为偶数,一般为2的幂次级,2,4,8,16...
k: 满足N=2^k(k>0),实质上k是N个采样数据可以分解为偶次幂和奇次幂的次数
real[]: inv=1时,存放N个采样数据的实部,inv=-1时,存放傅里叶变换的N个实部
imag[]: inv=1时,存放N个采样数据的虚部,inv=-1时,存放傅里叶变换的N个虚部
出口参数:
real[]: inv=1时,返回傅里叶变换的实部,inv=-1时,返回逆傅里叶变换的实部
imag[]: inv=1时,返回傅里叶变换的虚部,inv=-1时,返回逆傅里叶变换的虚部
*/
void FFT::dealFFT(double real[], double imag[], double dSp[], int N, int k, int inv)
{
int i, j, k1, k2, m, step, factor_step;
double temp_real, temp_imag, factor_real, factor_imag;
if (inv != 1 && inv != -1)
return;
//double *real = new double[N];
//double *imag = new double[N];
//倒序
j = 0;
for (i = 0; i < N; i++)
{
if (j>i)
{
temp_real = real[j];
real[j] = real[i];
real[i] = temp_real;
temp_imag = imag[j];
imag[j] = imag[i];
imag[i] = temp_imag;
}
m = N / 2;
while (j >= m&&m != 0)
{
j -= m;
m >>= 1;
}
j += m;
}
//蝶形运算
for (i = 0; i < k; i++)
{
step = 1 << (i + 1);
factor_step = N >> (i + 1); //旋转因数变化速度
//初始化旋转因子
factor_real = 1.0;
factor_imag = 0.0;
for (j = 0; j < step / 2; j++)
{
for (k1 = j; k1 < N; k1 += step)
{
k2 = k1 + step / 2; //蝶形运算的两个输入
/* temp_real = real[k1] + real[k2] * factor_real - imag[k2] * factor_imag;
temp_imag = imag[k1] + real[k2] * factor_imag + imag[k2] * factor_real;
real[k2] = real[k1] - (real[k2] * factor_real - imag[k2] * factor_imag);
imag[k2] = imag[k1] - (real[k2] * factor_imag + imag[k2] * factor_real);
real[k1] = temp_real;
imag[k1] = temp_imag;*/
//上面方法虽然直白,但效率太低,稍微改变结构如下:
temp_real = real[k2] * factor_real - imag[k2] * factor_imag;
temp_imag = real[k2] * factor_imag + imag[k2] * factor_real;
real[k2] = real[k1] - temp_real;
imag[k2] = imag[k1] - temp_imag;
real[k1] = real[k1] + temp_real;
imag[k1] = imag[k1] + temp_imag;
}
factor_real = inv*cos(-2 * PI*(j + 1)*factor_step / N);
factor_imag = inv*sin(-2 * PI*(j + 1)*factor_step / N);
}
}
if (inv == -1)
{
for (i = 0; i <= N - 1; i++)
{
real[i] = real[i] / N;
imag[i] = imag[i] / N;
}
}
for (i = 0; i<N;i++)
{
dSp[i] = sqrt(real[i] * real[i] + imag[i] * imag[i]);
}
}
一般好像需要进行下转换,即后半部分和前半部分置换,即1234变成3412.
void FFT::FFTShift(double dp[], int len)
{
for (int i = 0; i < len / 2; i++)
{
double tmp = dp[i];
dp[i] = dp[i + len / 2];
dp[i + len / 2] = tmp;
}
}
此时得到的应该是实部和虚部解出来的频谱图的Y轴电压值,一般频谱图Y轴为dB,因此需要进行转换
void FFT::getFFT(double dRe[], double dIm[], double dSp[], int len, int nBits, double dWorkingImpedance)
{
dealFFT(dRe, dIm, dSp, len, nBits, 1);
FFTShift(dSp,len); //此时得到的应该是实部和虚部解出来的频谱图的Y轴电压值,还需要转换
////dBW = 10lg(电压^2/阻抗);dBm =dBW+30,注意电压单位是V
for (int i = 0; i<len; i++)
{
dSp[i] = 10 * log10(dSp[i] * dSp[i] / dWorkingImpedance)+30;
}
}
getFFT()输出之后的dp才是要的频谱图Y轴值,频谱图X轴的坐标得到通过以下方式:
//X轴精确度,采样频率/数据个数 = 步长
m_DeltaX_S = m_dataPara.nSampleFrequency / nDataNumOfPage_S;
data_SX[i /2] = m_dataPara.nCenterFrequency + countm_DeltaX_S - m_dataPara.nWorkingBandWidth/2;//中心频率+当前点步长-带宽/2
在项目中,实际代码如下:
int count = 0;
for (int i = 0; i < nDataNumOfPage_S * 2; i++)
{
if (i % 2 == 0)
data_SQ[i / 2] = data_S[i] * m_DeltaY_S;
else
data_SI[i / 2] = data_S[i] * m_DeltaY_S;
if (i % 2 == 0)
{
count++;
data_SX[i / 2] = m_dataPara.nCenterFrequency + count*m_DeltaX_S - m_dataPara.nWorkingBandWidth/2;
}
}
m_dataPara.nWorkingImpedance = 50;
FFT fft;
int nBits = log10(nDataNumOfPage_S) / log10(2);//因为参数需要是2的N次方
fft.getFFT(data_SQ, data_SI, data_SS, nDataNumOfPage_S, nBits, m_dataPara.nWorkingImpedance); LoadData_S(data_SX, data_SS, nDataNumOfPage_S);
。。。
其他参考文章:
http://blog.sina.com.cn/s/blog_65d639d50101buo1.html
http://blog.csdn.net/hippig/article/details/8778753
http://www.makaidong.com/%E5%8D%9A%E5%AE%A2%E5%9B%AD%E6%8E%92%E8%A1%8C%E6%A6%9C/20151025/365773.html
原文链接: https://www.cnblogs.com/zwh0214/p/6297121.html
欢迎关注
微信关注下方公众号,第一时间获取干货硬货;公众号内回复【pdf】免费获取数百本计算机经典书籍
原创文章受到原创版权保护。转载请注明出处:https://www.ccppcoding.com/archives/247996
非原创文章文中已经注明原地址,如有侵权,联系删除
关注公众号【高性能架构探索】,第一时间获取最新文章
转载文章受原作者版权保护。转载请注明原作者出处!