引言
本文主要介绍对于音频信号进行谐波失真的Matlab仿真分析，到C语言在arm平台运行，并正确计算出失真度。
在开始本篇文章之前，结合之前文章的知识积累进行介绍。
首先是音频为什么会发生失真，失真的几个分类还有测量的方法
然后就是对音频信号中基波和谐波等概念进一步进行理解
有了这些准备，可以知道我们要想对信号进行失真度测试，那么就需要对信号进行FFT的操作，变换到频域进行分析。因为对于时域中的实际信号，对于我们而言是及其复杂的，无法进行直接的分析，傅里叶级数的出现，我们可以把所有的复杂信号都看成是正弦波组成。
通过傅里叶变换，我们就可以将时域中复杂的信息，变成我们关心的每个正弦信号的频率和幅度信息，从而有了对信号进行处理的可能。
其中关键点，一个是对信号进行时域到频域的转换，另一个更重要的是时域和频域中的频率和幅度究竟是上面样的关系。
从Matlab平台进行FFT到ARM平台C语言FFT频谱分析
FFT可以将一个信号的频谱提取出来，这在频谱分析方面也是经常用的。
虽然很多人都知道FFT是什么，可以用来做什么，怎么去做，但是却不知道FFT之后的结果是什么意思、如何决定要使用多少点来做FFT。所以对于越基础的概念，大家往往忘记的越厉害，其实答案早就都在课本中。
一个模拟信号，经过ADC采样之后，就变成了数字信号。采样定理告诉我们，采样频率要大于信号频率的两倍，这些我就不在此罗嗦了。采样得到的数字信号，就可以做FFT变换了。根据我们对FFT后对频率分辨率的要求，得到FFT的N个采样点，经过FFT之后，就可以得到N个点的FFT结果。为了方便进行FFT运算，通常N取2的整数次方。
Matlab谐波失真仿真分析
对于音频信号，最重要的指标便是谐波失真，对于失真的测量最常见的技术指标有谐波失真（HD）、总谐波失真（THD）、总谐波失真加噪声（Total Harmonic Distortion＋Noise，简称THD＋N）和互调失真（Intermodulation Distortion，简称IMD）。
我们这里就先来分析THD，也就是总谐波失真。
matlab THD函数解释
通过matlab的help，可以看到对于thd函数，我们需要用到的是
r = thd(x,fs,n) specifies the sample rate, fs, and the number of harmonics (including the fundamental) to use in the THD calculation.
详细解释如下
Specify Number of Harmonics
Open Live Script
Create a signal sampled at 1 kHz. The signal consists of a 100 Hz fundamental with amplitude 2 and three harmonics at 200, 300, and 400 Hz with amplitudes 0.01, 0.005, and 0.0025.
Set the number of harmonics to 3. This includes the fundamental. Accordingly, the power at 100, 200, and 300 Hz is used in the THD calculation.
t = 0:0.001:1-0.001;
x = 2*cos(2*pi*100*t)+0.01*cos(2*pi*200*t)+ ...
0.005*cos(2*pi*300*t)+0.0025*sin(2*pi*400*t);
r = thd(x,1000,3)
r = -45.0515
Specifying the number of harmonics equal to 3 ignores the power at 400 Hz in the THD calculation.
THD Matlab实现
既然是信号分析，有了理论依据，知道我们要干什么了，那么接下来就是通过matlab这种方便的工具进行快速的验证。
因为如果我们直接进行C语言代码的编写和应用，会导致我们会象个无头苍蝇，没有办法直观的知道哪一步应该是什么样的输出和结果，也没有办法快速的验证结果是不是正确，以及程序的理想程度。
matlab的仿真验证就非常的简单。主要步骤为
- 生成测试正弦信号
- THD函数验证
不多逼逼，先把代码贴上来
clc;close all;clear all;
M=8192;%fft采样点
Fs=48000; %采样频率，一秒多少个采样点
N=48000*2; %序列长度，总数据有多少个点，即对N个点FFT
f1=50;
f2=750;
f3=1500;
T=1/Fs; %单个点采样时间
t=(0:1:N-1)*T; %总时间T/Fs，每个点时间间隔t
y=2+6*sin(2*pi*f1*t-pi*30/180)+2*sin(2*pi*f2*t+pi*90/180)+2*sin(2*pi*f3*t+pi*60/180); %生成输入信号
% y=0.1*sin(2*pi*f1*t+pi/3)+0.2*sin(2*pi*f2*t-pi/6);
figure(1)
plot(t, y);
% 存储信号到wav文件
% filename = 'sine.wav';
% audiowrite(filename, y, Fs);
% Matlab FFT
Y = fft(y, M); %作FFT变换
A = abs(Y); %计算幅值
figure(2)
% plot(0:1:(N-1)/2, A(1:N/2));
plot(0:1:M-1, A(1:M));title('Matlab N=8192');xlabel('N');ylabel('Amplitude');%作图
grid on
% 总谐波失真TDH
NumHarmonics = 50; %谐波个数，保留的谐波个数
[thd_db, harmpow, harmfreq] = thd(y, Fs, NumHarmonics); %计算thd 单位为dB
percent_thd = 100 * (10^(thd_db / 20)); %转换为百分比
T = table(harmfreq, harmpow, 'VariableNames',{'Frequency','Power'});
figure(3);
thd(y, Fs, NumHarmonics); %画出谐波图
display( percent_thd); %打印最终总谐波失真（百分比）
display( thd_db);
fft的部分可以先不关心，先来分析THD的结果
上图中是THD的频率能量图，matlab程序中，选择的次谐波数量为50，其实不需要计算这么多的谐波个数，我们这个自己生成的信号，谐波数量为3，对于一般的音频信号谐波数量取5或者10就可以了，其他的一般可以忽略不计。
在命令行窗口还输出了一个百分比的结果，因为在音频领域，THD通常不用db表示而是通过百分比来表示。
我们通过maltab是来仿真程序是不是正确的，这里多说一嘴，那么我们又如何来验证我们的matlab计算的结果就是正确的呢？
最终还是绕不开通过我们最基本的公式来进行验算，总谐波失真是除基波外的全部的谐波幅值平方和开跟 比上 基波的幅值，谐波失真就是目标谐波幅值平方和开跟 比上 基波的幅值，对于matlab中的thd函数，其实更像是谐波失真，而不是总谐波失真，因此为了方便理解和计算，我们就直接使用已知的信号进行计算和验证，下文中的谐波失真和总谐波失真因为数据的已知性，所以都是一样的，请放心食用。
谐波计算
根据采样定理，当采样频率大于输入信号最高频率的2倍时，信号的频率可以无失真被恢复出来．信号的失真度定义为所有谐波能量之和比上基波能量的结果的平方根．当信号中的干扰远小于各次谐波的总能量时，失真度可由公式计算得出（此设计取5次谐波），即
基波功率P1可由频谱幅值U1求得．因此可写成
因为我们的测试信号就是我们自己生成的，所以其中采样率还有幅值等信息都是已知，对于我们本文要测试的信号y=2+6*sin(2*pi*f1*t-pi*30/180)+2*sin(2*pi*f2*t+pi*90/180)+2*sin(2*pi*f3*t+pi*60/180);，2是我们的直流分量，不在我们的谐波计算范围内，从基波开始到高次谐波 U1 = 6，U2 = 2， U3 = 2，直接带入公式就可以得到
THD = 根号下（U2平方 + U3平方） / U1 *100% = 根号下（2平方 + 2平方） / 6 = 2根号2 / 6 = 47.1405%
由此可以快速的验证matlab的thd计算得到的结果是正确的。
到了这里，其实也引出来了一个问题，对于一个未知的实际信号，我们不可能直接获得它各分量的幅值信息。但是我们有一个好帮手，就是频域分析，通过对采集到的信号进行FFT，就可以快速的得到采集信号的频域信息，从频域就可以简单方便的对信号的频率分量和幅值进行分析。
频域中得到基波以及各次谐波的幅值，带入我们的转换公式，可以得到对应的时域中的正弦分量中的幅值，从而就可以计算出该信号的失真度。
频域频率点对应时域正弦分量频率
假设某单一的正弦波（带有高阶谐波）y=2+6*sin(2*pi*f1*t-pi*30/180)+2*sin(2*pi*f2*t+pi*90/180)+2*sin(2*pi*f3*t+pi*60/180);，经过采样、FFT后得到它的频谱图如图所示．
假设FFT之后某点n用复数a+bi表示，那么这个复数的模就是An=根号（aa+bb），相位就是Pn=atan2(b,a)。根据以上的结果，就可以计算出n点（n≠1，且n<=N/2）对应的信号的表达式为：
An / (N/ 2)* cos(2* pi* Fn* t+Pn)，
即 2* An/N* cos(2* pi* Fn* t+ Pn)。
我们现在已经很清晰的知道时域中要被fft的信号
频率=0 ，直流分量 ；幅值=2v；相移0；
频率=50hz，基波；幅值=6v；相移-30；
频率=750hz，2次谐波；幅值=2v；相移+90；
频率=1500hz，3次谐波；幅值=2v；相移+60；
现在就是通过频域来进行反向推算和验证。
首先是对频率进行对号入座
信号的采样频率为48khz，fft的采样点为8192，所以，fft的频率分辨率为 ：信号采样频率 / FFT采样点 = 48k / 8192 = 5.859375，约等于6，也就是fft的频率误差可能在6个频率点左右
对于频谱图中的0点，即fft的第1个点，对应的就是时域中直流分量，也就是0次谐波，0 * 5.859375 = 0，所以时域中频率为0
对于频谱图中的9点，即fft的第10个点，对应的就是信号中的主频（基频），一次谐波，9 * 5.859375 = 52.734375，所以对应时域中频率为50hz的正弦分量
对于频谱图中的128点，即fft的第129个点，对应的就是时域中的f2，高频谐波，二次谐波，128 * 5.859375 = 750，所以对应时域中频率为750hz的正弦分量
对于频谱图中的258点，即fft的第259个点，对应的就是时域中的f3，高频谐波，三次谐波，256 * 5.859375 = 1500，所以对应时域中频率为1500hz的正弦分量
/ ******************************** /
2021/11/27 更正，图中2次谐波 3次谐波表述错误，而是15次谐波和100次谐波，这种高次谐波通常不会被自然的非线性谐波失真产生，也不会有这么大的幅值，这里仅仅是为了说明时域和频域的关系。
/ ******************************** /
频域谐波幅值对应时域正弦分量幅值
频谱图中的横坐标解释完毕，现在开始对频谱中的幅值进行分析，频域中幅值和时域中幅值连接起来的变量是FFT的采样点
直流分量
频域幅值 = 时域幅值 * FFT采样点
正弦分量
频域幅值 = 时域幅值 * FFT采样点/2
所以
频谱图中0点，幅值为18047.35，则时域中对应频率为0的正弦分量，即直流分量的幅值为 18047.35 / 8192 = 2.20304；
频谱图中9点，幅值为16379.8909，则时域中对应频率为50hz的正弦分量，即基波的时域幅值为 16379.8909 / （8192/2） =3.9989968017578125；
频谱图中128点，幅值为8182.281，则时域中对应频率为750hz的正弦分量，即2次谐波的时域正弦分量的幅值为 8182.281 / （8192/2）= 1.997627197265625；
频谱图中256点，幅值为8175.976，则时域中对应频率为1500hz的正弦分量，即3次谐波的时域正弦分量的幅值为 8175.976 / （8192/2） = 1.996087890625；
依然是因为FFT采样点我们选择的是8192，但是信号的采样频率为48000，所以导致分辨率不足，出现了FFT运算中常见的频率泄露现象。我们的基波幅值应该是6，而不是计算得到的3.9989968017578125。所以，我们要对FFT的分辨率进行提高。
我们把FFT的采样点变成 8192*4 = 32768
即分辨率为 Fs / 32768 = 1.46
clc;close all;clear;
M=8192*4;%fft采样点
Fs=48000; %采样频率，一秒多少个采样点
N=48000*2; %序列长度，总数据有多少个点，即对N个点FFT
f1=50;
f2=750;
f3=1500;
T=1/Fs; %单个点采样时间
t=(0:1:N-1)*T; %总时间T/Fs，每个点时间间隔t
y=2+6*sin(2*pi*f1*t-pi*30/180)+2*sin(2*pi*f2*t+pi*90/180)+2*sin(2*pi*f3*t+pi*60/180); %生成输入信号
% y=6*sin(2*pi*f1*t+pi/3);
figure(1)
plot(t, y);
% 存储信号到wav文件
% filename = 'sine.wav';
% audiowrite(filename, y, Fs);
% Matlab FFT
Y = fft(y, M); %作FFT变换
A = abs(Y); %计算幅值
figure(2)
% plot(0:1:(N-1)/2, A(1:N/2));
plot(0:1:M-1, A(1:M));title('Matlab N=2^15');xlabel('N');ylabel('Amplitude');%作图
grid on
/ ******************************** /
2021/11/27 更正，图中2次谐波 3次谐波表述错误，而是15次谐波和100次谐波，这种高次谐波通常不会被自然的非线性谐波失真产生，也不会有这么大的幅值，这里仅仅是为了说明时域和频域的关系。
/ ******************************** /
现在我们再次进行计算幅值
频谱图中0点，幅值为65455.6，则时域中对应频率为0的正弦分量，即直流分量的幅值为 65455.6 / 32768= 1.99754638671875；
频谱图中9点，幅值为95271.9，则时域中对应频率为50hz的正弦分量，即基波的时域幅值为 95271.9 / （32768/2） = 5.814935302734375；
频谱图中128点，幅值为32766，则时域中对应频率为750hz的正弦分量，即2次谐波的时域正弦分量的幅值为 32766 / （32768/2）= 1.9998779296875；
频谱图中256点，幅值为32753.6，则时域中对应频率为1500hz的正弦分量，即3次谐波的时域正弦分量的幅值为 32753.6 / （32768/2） = 1.99912109375；
导致最后的结果与我们的理论值有一定的差别
其实到了这里，我们fft的目的也就达成了，就是为了通过频域可以更加直观的分析出时域中各正弦分量的幅值，并且对其带入总谐波失真的公式进行计算。
那么，我们就先带入我们FFT后反推得到的幅值运算，看看与我们的理论值可以差多少。
t h d = U 2 2 + U 3 2 U 1 = 1.9998779296875 2 2 + 1.99912109375 6 2 5.814935302734375 thd=\frac{\sqrt{{{U}_{2}}^{2}+{{U}_{3}}^{2}}}{{{U}_{1}}}=\frac{\sqrt{\text{1}\text{.9998779296875}{{\text{2}}^{2}}+\text{1}\text{.99912109375}{{\text{6}}^{2}}}}{\text{5}\text{.814935302734375}} thd=U1U22+U32=5.8149353027343751.999877929687522+1.9991210937562
计算结果为0.48628561297026372798437583004506
我们的理论值为0.471405，误差在5%以内，如果提高FFT的采样点，结果将接近我们的理论值
发现知乎一大佬对于这块讲解比较详细，如果还是有疑问的可以跳转一下链接
https://zhuanlan.zhihu.com/p/85863024
C语言谐波失真应用
有了matlab的仿真，帮助我们理解了为什么要进行FFT，以及我们变换到频域后要做什么。有了这些作为操纵的依据，那么我们可以进行C代码的功能逻辑编写了。
首先梳理一下步骤
- 对采样信号的N点进行FFT
- 频域中分量分析，换算时域分量的幅值
- 计算THD
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define PI 3.1415926535
void kfft(double pr[], double pi[], int n, int k, double fr[], double fi[])
{
int it, m, is, i, j, nv, l0;
double p, q, s, vr, vi, poddr, poddi;
for (it = 0; it <= n-1; it++) { //将pr[0]和pi[0]循环赋值给fr[]和fi[]
m = it;
is = 0;
for(i = 0; i <= k-1; i++) {
j = m / 2;
is = 2 * is + (m - 2 * j);
m = j;
}
fr[it] = pr[is];
fi[it] = pi[is];
}
pr[0] = 1.0;
pi[0] = 0.0;
p = 6.283185306 / (1.0 * n);
pr[1] = cos(p); //将w=e^-j2pi/n用欧拉公式表示
pi[1] = -sin(p);
for(i = 2; i <= n-1; i++) { //计算pr[]
p = pr[i - 1] * pr[1];
q = pi[i - 1] * pi[1];
s = (pr[i - 1] + pi[i - 1]) * (pr[1] + pi[1]);
pr[i] = p - q;
pi[i] = s - p - q;
}
for(it = 0; it <= n-2; it = it + 2) {
vr = fr[it];
vi = fi[it];
fr[it] = vr + fr[it + 1];
fi[it] = vi + fi[it + 1];
fr[it + 1] = vr - fr[it + 1];
fi[it + 1] = vi - fi[it + 1];
}
m = n / 2;
nv = 2;
for(l0 = k - 2; l0 >= 0; l0--) { //蝴蝶操作
m = m / 2;
nv = 2 * nv;
for(it = 0; it <= (m - 1) * nv; it = it + nv)
for (j = 0; j <= (nv / 2) - 1; j++) {
p = pr[m * j] * fr[it +j + nv / 2];
q = pi[m * j] * fi[it +j + nv / 2];
s = pr[m * j] + pi[m * j];
s = s * (fr[it + j + nv / 2] + fi[it + j + nv / 2]);
poddr = p - q;
poddi = s - p - q;
fr[it + j + nv / 2] = fr[it + j] - poddr;
fi[it + j + nv / 2] = fi[it + j] - poddi;
fr[it + j] = fr[it + j] + poddr;
fi[it + j] = fi[it + j] + poddi;
}
}
for(i = 0; i <= n - 1; i++) {
pr[i] = sqrt(fr[i] * fr[i] + fi[i] * fi[i]); //幅值计算
}
return;
}
int main()
{
int i,j;
int f1 = 50;
int f2 = 750;
int f3 = 1500;
int N = 48000*2;
int Fs = 48000;
double T = 1.0 / Fs;
double pr[N], pi[N], fr[N], fi[N], t[N];
FILE *fp_fftdata, *fp_sine; //文件指针
fp_fftdata=fopen("fftdata.txt","w"); //fopen打开文件，这个文件可以是当前不存在的。“w”以写入的形式打开，“r”以读的形式打开
if(fp_fftdata==NULL) { //判断如果文件指针为空
printf("File cannot open! " );
exit(0); //在以0的形式退出，必须在文件开头有#include <stdlib.h>,stdlib 头文件即standard library标准库头文件
}
fp_sine=fopen("sinedata.txt","w"); //fopen打开文件，这个文件可以是当前不存在的。“w”以写入的形式打开，“r”以读的形式打开
if(fp_sine==NULL) { //判断如果文件指针为空
printf("File cannot open! " );
exit(0); //在以0的形式退出，必须在文件开头有#include <stdlib.h>,stdlib 头文件即standard library标准库头文件
}
for (i=0; i<N; i++) { //生成输入信号
t[i] = i*T;
pr[i] = 2+6*sin(2*PI*f1*t[i]-PI/6)+2*sin(2*PI*f2*t[i]+PI/2)+2*sin(2*PI*f3*t[i]+PI/3)+2*sin(2*PI*1600*t[i]+PI/3)+2*sin(2*PI*1800*t[i]+PI/3);
pi[i] = 0.0;
fprintf(fp_sine,"%f\n", pr[i]);
// printf("%d\t%lf\n",i,pr[i]); //输出结果
}
int M = 8192*4;
kfft(pr,pi,M,15,fr,fi); //调用FFT函数 fft取样点M = 2^k
for (i=0; i<M; i++) {
fprintf(fp_fftdata,"%f\n", pr[i]); //写入指针fp_fftdata，写入的东西就是刚才的用户输入的d,注意这里的fp_fftdata和d没有引号
// printf("%d\t%f\n",i,pr[i]); //输出结果
}
double direct_frequency = 0, direct_amplitude = 0;
if(pr[0] > pr[1]) {
direct_frequency = pr[0] * 0;
direct_amplitude = pr[0] / M;
}
int z = 0, flag = 0;
double fft_data_n[100] = {0},fft_data_a[100] = {0};
double molecule = 0.0,denominator = 0.0,thd = 0.0;
for(i = 2; i < M/2 - 1; i++) {
if(pr[i] >= pr[i-1] && pr[i] >= pr[i+1]) {
fft_data_n[z] = i;
fft_data_a[z] = pr[i];
z++;
}
}
denominator = fft_data_a[0] / (M / 2.0);
for(i = 1; i<z; i++) {
molecule = pow(fft_data_a[i] / (M / 2.0), 2) + molecule;
}
molecule = sqrt(molecule);
thd = molecule / denominator;
printf("thd is:%f\n", thd);
//关闭文件
fclose(fp_fftdata);
fclose(fp_sine);
return 0;
}
通过maltab仿真测试得到的THD结果为
y=2+6*sin(2*pi*f1*t-pi*30/180)+2*sin(2*pi*f2*t+pi*90/180)+2*sin(2*pi*f3*t+pi*60/180); %生成输入信号
% 总谐波失真TDH
NumHarmonics = 30; %谐波个数，保留的谐波个数
[thd_db, harmpow, harmfreq] = thd(y, Fs, NumHarmonics); %计算thd 单位为dB
percent_thd = 100 * (10^(thd_db / 20)); %转换为百分比
T = table(harmfreq, harmpow, 'VariableNames',{'Frequency','Power'});
figure;
thd(y, Fs, NumHarmonics); %画出谐波图
display( percent_thd ); %打印最终总谐波失真（百分比）
display( thd_db );
看似相差不大，但是对于失真度测试来讲1%的误差还是不能接受的，并且实际信号也不会这么完美，而是更多的16位、24位位深的定点数音频信号。
因此，还是要继续对算法程序进行改进。
8306
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
