C语言总谐波失真（THD）实现，从理论到应用分析改进详解
Posted “逛丢一只鞋”
tags:
篇首语：本文由小常识网(cha138.com)小编为大家整理，主要介绍了C语言总谐波失真（THD）实现，从理论到应用分析改进详解相关的知识，希望对你有一定的参考价值。
引言
从Matlab总谐波失真（THD）仿真到C语言总谐波失真（THD）应用
对于如何实现THD，上篇文章中已经叙述的比较清晰，但是，正如结尾中表述，实际计算数据与理论数据差距过大，无法应用在实际的系统中。
测试信号生成
为了更好的分析THD，依然是生成一个测试正弦信号
通过matlab生成，然后再去使用au软件进行处理
M=8192*8;%fft采样点
Fs=48000; %采样频率，一秒多少个采样点
N=48000*4; %序列长度，总数据有多少个点
f1=100;
f2=200;
f3=300;
T=1/Fs; %单个点采样时间
t=(0:1:N-1)*T; %总时间T/Fs，每个点时间间隔t
y=0.1*sin(2*pi*f1*t-pi/6)+0.0001*sin(2*pi*f2*t+pi/2)+0.0001*sin(2*pi*f3*t+pi/3); %生成输入信号
% % y=6*sin(2*pi*f1*t+pi/3);
% figure(1)
% plot(t, y);
% y = awgn(y, 0, 'measured');
% 存储信号到wav文件
filename = 'sine_24bit_48k.wav';
audiowrite(filename, y, Fs, 'BitsPerSample', 24);
[y, Fs] = audioread('sine_24bit_48k.wav'); %读取wav文件
figure
plot(t, y);
Y = fft(y, M); %作FFT变换
A = abs(Y); %计算幅值
figure
% plot(0:1:(N-1)/2, A(1:N/2));
plot(0:1:M-1, A(1:M));title('Matlab N=131,072');xlabel('N');ylabel('Amplitude');%作图
最终生成一个24位，48khz的wav文件
wav音频文件解析读取 定点转浮点分析 幅值提取（C语言实现）
正如此篇文章，对于wav文件解析以及提取就不再赘述
THD分析
采样点选择
对于THD，原理上没有什么可深挖的，就是谐波能量与基波能量的比值，公式变形后就是谐波分量的幅值平方和相加再开根，然后比上基波的幅值
对与理想的信号源以及FFT后，应该出现的是如上图所示的频谱图，所以在我们之前的C语言代码中，简单粗暴的只要频谱图中这个信号是高于两侧信号，就认为这个是谐波信号。
但是，对于这么一个24位的wav音频文件，我们得到的实际fft后信号没有这么的完美。
罪魁祸首就是采样点导致的频率分辨率的问题
我们这段信号采用频率为48000Hz，FFT的采样点为 8192*8 = 65536，按理来说，频率分辨率为 48000 / 65536 = 0.732421875
所以，对于正弦信号来讲基波的理论频域频率点应该是 100 / 0.732421875 = 136.53333333333333333333333333333
理论频域中的幅值应该是 0.1 * 65536 / 2 = 32768
但是实际基波的值为 频率点 137，频率幅值 2220.34
频率点值还好，差别不大，可以忽略误差，但是对于频率的幅值来讲，简直就是灾难，如果使用实际中基波频域的幅值去换算时域的幅值，得到的是 2220.34 / 65536 * 2 = 0.067759399 ，与理论值的0.1幅值相比相差甚远。
除了这个问题，在采样点为65536的程序中还存在一个BUG
例如在频谱图中的点273处
我们可以看到不同于我们一般见到的谐波信号是一个脉冲，这里是一个凹陷
但是，根据谐波是基波频率的整数倍的原理，我们知道在273频点上就是我们的二次谐波
这个其实也是采样频率0.7带来的偏差，虽然比采样频率1更加精准了，但是也会导致一些点的遗漏
改进后
所以，针对这个问题，解决方法就在FFT采样点的选择上
在matlab仿真的过程中，采样点选择65536，按理来说处理一个采样频率为48000采样率的信号绰绰有余，但是现实就是无情的打脸（基础不牢，受制于C语言中基2的FFT算法）
因此，增加FFT的采样点，这一次增加到 8192*16 = 131,072
现在频率的分辨率来到了 48000 / 131072= 0.3662109375
这时候我们再来看频谱图
现在对这个采样点又增加了一倍的正弦信号分析
对于正弦信号来讲基波的理论频域频率点应该是 100 / 0.3662109375 = 273.0666666666666
理论频域中的幅值应该是 0.1 * 131072/ 2 = 6,553.6
但是实际基波的值为 频率点 273，频率幅值 6505.15，根据实际基波逆推得到的时域频率为 99.9755859375，时域幅值为 0.099260711669921875
这一次，我们终于正儿八经的接近了正确的答案，得到了理想的值
所以带着matla的仿真结果，就可以去修改C语言中FFT的采样点，以及判断C语言中各阶段结果与Matlab结果的偏差
代码实现
有了充足的理论支撑，我们可以对C语言进行大刀阔斧的动手了
需要注意的有以下几点
- 基波的判断，要先把所有的频域分量找出来，然后找到幅值最大的那个基波，一个是得到频域的幅值，另一个就是得到频域的频率点
- 谐波的判断，首先谐波要在基波时域频率的整数倍上，其次谐波除了要比两边的频率点幅值大，还要保证不是误差引起的只是大了一丢丢的意外
- 噪声的判断，这里直接把噪声也考虑了进来，对于噪声来说，也要确保是真正的噪声，最简单的方法就是比较与周围频率的幅值
#define PI acos(-1)
#define N 48000*4 // 输入数据长度
#define Fs 48000
#define M 8192*8 // FFT 采样点 2的n次幂
#define Mn 16 // FFT 采样点次幂数 2^Mn
void THD_TEST_Pthread(void)
while(1)
if(fft_flag == 1)
int i,j;
double T = 1.0 / Fs;
double pr[N];
double pi[N], fr[N], fi[N], t[N];
int z = 0;
int fft_data_n[256], denominator_n = 0;
double fft_data_a[256];
double molecule = 0.0, denominator = 0.0, thd = 0.0;
double direct_frequency = 0.0, direct_amplitude = 0.0;
int temp = 0;
char thd_data[4];
int ret;
double *wavdata;
FILE *fp_fftdata, *fp_sine; //文件指针
fp_fftdata=fopen("fftdata.txt","w");
if(fp_fftdata==NULL)
printf("File cannot open! " );
exit(0);
fp_sine=fopen("sinedata.txt","w");
if(fp_sine==NULL)
printf("File cannot open! " );
exit(0);
// 初始化要使用的数组
memset(pr, 0.0, sizeof(pr));
memset(pi, 0.0, sizeof(pi));
memset(fr, 0.0, sizeof(fr));
memset(fi, 0.0, sizeof(fi));
memset(t, 0.0, sizeof(t));
memset(fft_data_n, 0.0, sizeof(fft_data_n));
memset(fft_data_a, 0.0, sizeof(fft_data_a));
memset(thd_data, 0, sizeof(thd_data));
fft_flag = 0;
// int f1 = 50;
// int f2 = 750;
// int f3 = 1500;
// for (i = 0; i < N; i++) //生成输入信号
// t[i] = i * T;
// // pr[i] = 2 + 6 * sin(2 * PI * f1 * t[i] - PI / 6) + 0.02 * sin(2 * PI * f2 * t[i] + PI / 2) + 0.01 * sin(2 * PI * f3 * t[i] + PI / 3);
// pr[i] = 0.1*sin(2 * PI * f1 * t[i] - PI / 6) + 0.0001 * sin(2 * PI * f2 * t[i] + PI / 2) + 0.0001 * sin(2 * PI * f3 * t[i] + PI / 3);
// pi[i] = 0.0;
// fprintf(fp_sine,"%f\\n", pr[i]);
// // printf("%d\\t%f\\n",i,pr[i]); //输出结果
//
wavdata = (double *) calloc ( N, sizeof(double) );
ret = parse_wave_file_double("sine_24bit_48k.wav", wavdata);
if(ret != 0)
ERRO("parse_wave_file Error \\n");
break;
for(i = 0; i < N; i++)
pr[i] = wavdata[i];
fprintf(fp_sine, "%.15f\\n", pr[i]);
free(wavdata);
// FFT
kfft(pr, pi, M, Mn, fr, fi); //调用FFT函数 fft取样点M = 2^k
for (i = 0; i < M; i++)
fprintf(fp_fftdata, "%.15f\\n", pr[i]); //写入指针fp_fftdata
// printf("%d\\t%f\\n",i,pr[i]); //输出结果
// 直流分量
if(pr[0] > pr[1])
direct_frequency = pr[0] * 0.0;
direct_amplitude = pr[0] / (M * 1.0);
z = 0;
molecule = 0.0;
denominator = 0.0;
denominator_n = 0;
thd = 0.0;
temp = 0;
// 谐波分量
for(i = 1; i < (M/2 - 1); i++)
if( (pr[i] > pr[i - 1]) && (pr[i] > pr[i + 1]) )
if(pr[i] / pr[i - 1] >= 2)
fft_data_n[z] = i;
fft_data_a[z] = pr[i];
z++;
if(z > 255)
z = 0;
break;
int fft_len = (double) sizeof(fft_data_a) / sizeof(*fft_data_a);
bubble_sort(fft_data_a, fft_data_n, fft_len);
denominator = fft_data_a[0] / (M / 2.0); // 基波分量 时域幅值
double Fs_Resolution = 0.0; // 频率分辨率
Fs_Resolution = Fs / (M*1.0);
denominator_n = (int)(fft_data_n[0] * Fs_Resolution + 0.5); // 基波分量 时域频率
int count_n = 0;
double fft_data[256];
memset(fft_data, 0, sizeof(fft_data));
for(i = 1; i < 255; i++)
temp = (int)(fft_data_n[i] * Fs_Resolution + 0.5);
if( temp % denominator_n == 0 )
fft_data[count_n] = fft_data_a[i];
count_n++;
if(count_n > 255)
count_n = 0;
break;
for(i = 0; i < count_n; i++)
molecule = pow(fft_data[i] / (M / 2.0), 2) + molecule;
molecule = sqrt(molecule); // 各次谐波平方和开根
thd = molecule / denominator;
printf("THD is : %f(percentage)\\n", thd);
temp = thd * 100000000;
//int数写入char数组
// char *thd_char = (char *)&temp;
for (i = 0; i < 4; i++)
thd_data[i] = (temp >> (8 * i)) & 0xff;
tcp_send_qt(0x0a, thd_data, 4);
//关闭文件
fclose(fp_fftdata);
fclose(fp_sine);
以上是关于C语言总谐波失真（THD）实现，从理论到应用分析改进详解的主要内容，如果未能解决你的问题，请参考以下文章