首页 > 其他分享 >音频信号频谱分析

音频信号频谱分析

时间:2022-10-15 21:35:26浏览次数:55  
标签:频谱 real struct 音频 complex 3.2 信号 imag

1 音频信号频谱分析简介

1.1 设计背景

在当今数字化时代的背景下,由于DSP具有处理速度快,功耗低,性能好,存储容量大等特点,因此已被广泛应用于通信、计算机、消费类电子产品等多个领域中。

而对于音频信号的分析来说,单纯只对其时域进行分析,有时候往往会不尽如人意,而如果换一种角度,从频域上对其进行分析,有时可能会有惊人的效果。如在一段嘈杂的音频中对其频域进行分析,便有可能从中确定声音的来源,再通过对其进行相应的滤波处理便能够很好地提取到对我们有用的信息。

因此通过DSP对音频信号进行频谱分析将会具有很好的实际意义和应用前景。

 

1.2 设计方法

本次课程设计将采用AIC23音频处理芯片完成对音频信号的采集,同时通过TMS320F28335对采集进来的音频信号进行FFT处理并将结果显示于CCS中从而完成对音频信号的频谱分析。

 

1.3 设计要求

通过A/D采集音频信号(fmax<10kHz),编写DSP的FFT处理程序(自定频谱分辨力),获得幅频特性之后,在点阵液晶中大致显示出幅频图。并在液晶中用文字显示频率幅值前三的频率值。

1、DSP与A/D接口电路的原理图绘制;

2、DSP控制A/D的程序编写与调试;

3、编写DSP的FFT处理程序;

4、控制点阵液晶,实现绘图功能,将幅频图显示出来

5、按要求编写课程设计报告书,正确、完整的阐述设计和实验结果。

6、在报告中绘制程序的流程图,并文字说明。

 

2 硬件设计

2.1 TMS320F28335最小系统

TMS320F28335最小系统是指能使DSP正常工作的最简电路,如图2.1所示,其中大体包括以下四个部分:

  1. 主芯片;
  2. 晶振电路模块:外接30MHz的晶振,内部可通过锁相环PLL将其倍频至150MHz作为系统的时钟频率;
  3. 复位电路模块:复位时低电平有效。上电后,由于电容两端电压不能突变的缘故,所以XRS端会在很短的一段时间内保持低电平,再被拉到高电平。复位后,DSP会先响应复位中断信号进入到中断向量表中,再通过中断向量表的指引执行BootLoader程序,完成对DSP的初始化工作;
  4. JTAG电路模块:用于下载程序和对程序进行仿真。

图2.1 TMS320F28335最小系统

 

2.2 音频编解码模块

2.2.1 TLV320AIC23芯片介绍

TLV320AIC23是一个高性能的多媒体数字语音编解码器,它的内部ADC和DAC转换模块带有完整的数字滤波器。详细指标如下:

①在 ADC 采用 96kHz 采样率时噪音 90-dB;

②在 DAC 采用 96kHz 采样率时噪音 100-dB;

③ 1.42V – 3.6V 核心数字电压:兼容 TI F28x DSP 内核电压;

④ 2.7V – 3.6V 缓冲器和模拟:兼容 TI F28x DSP 内核电压;

⑤支持 8kHz – 96kHz 的采样频率;

⑥软件控制通过 TI IIC 接口;

⑦音频数据输入输出通过 TI McBSP 接口。

 

2.2.2 TLV320AIC23硬件连接

如图2.2.2所示,LINE_IN2为音频输出口,可接耳机或音响播放音频;LINE_IN3为语音输入口,可接话筒。同时F28335通过IIC接口和McBSP接口对TLV320AIC23进行控制。

图2.2.2 TLV320AIC32硬件连接

 

3 软件设计

3.1 音频编解码模块初始化

如表3.1.1所示,根据芯片的Register Map,可通过IIC协议修改芯片内对应寄存器的值。

表3.1.1 TLV320AIC23芯片Register Map

ADDRESS

REGISTER

0000000

Left line input channel volume control

0000001

Right line input channel volume control

0000010

Left channel headphone volume control

0000011

Right channel headphone volume control

0000100

Analog audio path control

0000101

Digital audio path control

0000110

Power down control

0000111

Digital audio interface format

0001000

Sample rate control

0001001

Digital interface activation

0001111

Reset register

如表3.1.2所示,再根据芯片手册所提供的对"Sample rate control"的设定,即可将采样速率设置为8kHz。至此,音频编解码模块的初始化工作便已基本结束。

表3.1.2 TLV320AIC23芯片Sample rate control

SAMPLING RATE

FILTER TYPE

SAMPLING-RATE CONTROL SETTINGS

ADC
(kHz)

DAC
(kHz)

SR3

SR2

SR1

SR0

BOSR

96

96

3

0

1

1

1

0

88.2

88.2

2

1

1

1

1

1

48

48

0

0

0

0

0

0

44.1

44.1

1

1

0

0

0

1

32

32

0

0

1

1

0

0

8.021

8.021

1

1

0

1

1

1

8

8

0

0

0

1

1

1

48

8

0

0

0

0

1

0

44.1

8.021

1

1

0

0

1

1

8

48

0

0

0

1

0

0

8.021

44.1

1

1

0

1

0

1

 

3.2 音频编解码模块中断程序

由于本次FFT采用的点数是256,点数较多,运算量较大,因此不可能每次进入中断后都对音频进行FFT的处理,因此为了解决这个问题,我通过将每次采集到的音频都用FIFO的方法放入一个数组中储存起来(也即是队列),再判断新增音频的数量(也即是进入中断的次数),如果新增音频的数量等于FFT点数则对整个数组进行FFT处理并将记录新增音频数的变量清0;如果不是,则退出中断,如图3.2所示。这样便能避免每新增一个音频点便进行一次FFT计算这样低效率的冗余操作了。

图3.2 音频编解码模块中断程序流程图

下面对FFT程序进行介绍。注:FFT程序是先在Visual Studio中利用图形化界面调试好,再移植到CCS中的。

 

3.2.1 基本复数运算的定义

由于FFT程序涉及到大量的复数运算,因此要想进行FFT,首先最基本的便是定义基本的复数运算,如程序3.2.1所示,这里我们分别定义了复数的加法运算、减法运算和乘法运算。

 1 struct complex c_add(struct complex a, struct complex b)//复数加法
 2 {
 3     struct complex c;
 4     c.imag = a.imag + b.imag;
 5     c.real = a.real + b.real;
 6     return c;
 7 }
 8 struct complex c_sub(struct complex a, struct complex b)//复数减法
 9 {
10     struct complex c;
11     c.imag = a.imag - b.imag;
12     c.real = a.real - b.real;
13     return c;
14 }
15 struct complex c_mul(struct complex a, struct complex b)//复数乘法
16 {
17     struct complex c;
18     c.real = a.real*b.real - a.imag*b.imag;
19     c.imag = a.real*b.imag + a.imag*b.real;
20     return c;
21 }

程序3.2.1 基本复数运算的定义程序

 

3.2.2 旋转因子的定义

旋转因子的数学表达式如下,根据该表达式便可写出相应函数,如程序3.2.2所示。

    $W_{N}^{kn}=e^{j\left( -\frac{2\pi}{N} \right) kn}$

 

 1 struct complex get_WN_kn(double n, double N, double k)//获得旋转因子
 2 {
 3     int i = 0;
 4     struct complex ex, W;
 5     ex.real = 1.0*cos(2 * PI / N);
 6     ex.imag = -1.0*sin(2 * PI / N);
 7     W.real = 1; W.imag = 0;
 8     for (i = 0; i < k*n; i++)
 9         W = c_mul(ex, W);
10     return W;
11 }

程序3.2.2 旋转因子的定义程序 

 

3.2.3 码位颠倒函数

由于在FFT运算流图中,第一步便是把输入序列以倒序的形式输入,因此如程序3.2.3所示,该函数的大致思路是将原码通过右移的方法一位一位地送到倒序的码位上,最终实现码位颠倒,如图3.2.3所示。

图3.2.3 码位颠倒示意图

 1 void change(struct complex x[])
 2 {
 3     int i = 0, j = 0, k = 0;
 4     double t;
 5     struct complex temp;
 6     for (i = 0; i<fft_N; i++)
 7     {
 8         k = i; j = 0;
 9         t = (log(fft_N) / log(2));
10         while ((t--)>0)    //利用按位与以及循环实现码位颠倒
11         {
12             j = j << 1;
13             j |= (k & 1);
14             k = k >> 1;
15         }
16         if (j>i)    //将x(n)的码位互换
17         {
18             temp = x[i];
19             x[i] = x[j];
20             x[j] = temp;
21         }
22     }
23 }

程序3.2.3 码位颠倒函数程序

 

 

3.2.4 FFT之蝶形运算

首先调用码位颠倒函数,再根据蝶形运算的方法和流程进行FFT的最终运算,如图3.2.4和程序3.2.4所示。

图3.2.4 蝶形运算示意图

 1 void fft(struct complex x[], int N)
 2 {
 3     change(x);
 4     struct complex W, up, down, product;
 5     int i, l, j, k;
 6     for (i = 0; i<log(fft_N) / log(2); i++)   /*一级蝶形运算*/
 7     {
 8         l = 1 << i;
 9         for (j = 0; j<fft_N; j += 2 * l)      /*一组蝶形运算*/
10         {
11             for (k = 0; k<l; k++)             /*一个蝶形运算*/
12             {
13                 W = get_WN_kn(1, fft_N, k * fft_N / (2 * l));
14                 product = c_mul(x[j + k + l], W);
15                 up = c_add(x[j + k], product);
16                 down = c_sub(x[j + k], product);
17                 x[j + k] = up;
18                 x[j + k + l] = down;
19             }
20         }
21     }
22 }

程序3.2.4 FFT之蝶形运算程序

 

3.2.5 IFFT函数

该函数在本次课设中其实并没有起到实际的作用,但是在调试过程中却能用于检验FFT函数的正确性,结果如图3.2.5所示。

图3.2.5 FFT&IFFT函数在VS上的运行结果

 1 void ifft(struct complex x[], int N)
 2 {
 3     int i;
 4     for (i = 0; i<fft_N; i++)
 5         x[i].imag = -x[i].imag;
 6     fft(x, fft_N);
 7     for (i = 0; i<fft_N; i++)
 8     {
 9         x[i].real = x[i].real*1.0 / fft_N;
10         x[i].imag = -x[i].imag*1.0 / fft_N;
11     }
12 }

程序3.2.5.2 IFFT函数程序

 

4 运行结果

通过XDS100v2仿真器运行程序,并使用CCS中的画图工具,对经过FFT处理后的数组数据进行绘制便可得到以下结果。

4.1 "1000Hz+2500Hz"正弦波运行结果

如图4.1所示,该段音频是通过Matlab编写函数来产生的,Matlab函数请见附录A。

图4.1 "1000Hz+2500Hz"正弦波运行结果

 

4.2 "1500Hz+2000Hz"正弦波运行结果

如图4.2所示,该段音频是通过Matlab编写函数来产生的,Matlab函数请见附录A。

图4.2 "1500Hz+2000Hz"正弦波运行结果

 

4.3 "1700Hz+1800Hz"正弦波运行结果

如图4.3所示,该段音频是通过Matlab编写函数来产生的,Matlab函数请见附录A。

图4.3 "1700Hz+1800Hz"正弦波运行结果

 

4.4 "光辉岁月"音乐运行结果

如图4.4所示,光辉岁月这首歌的频率大多集中在342Hz这个频率左右。

图4.4 "光辉岁月"音乐运行结果

 

4.5 "Opera 2"音乐运行结果

如图4.5所示,Opera 2这首歌的最高频率可以达到2090Hz。

图4.5 "Opera 2"音乐运行结果

 

4.6 "相爱很难"音乐运行结果

如图4.6所示,相爱很难这首歌包含了两个主要频率,男音部分的频率大概是228Hz左右,女音部分的频率大概是532Hz左右。

图4.6 "相爱很难"音乐运行结果

 

附录A Matlab源程序代码

f1=600;%正弦频率
f2=200;%正弦频率
fs=22050;%采样频率
T=5;%总时长
x=1:T*fs;
y1=sin(f1*2*pi*x/fs);
y2=sin(f2*2*pi*x/fs);
y=y1;
sound(y,fs)

 

附录B DSP源程序代码

  1 #include "math.h"
  2 #define PI 3.1415926
  3 #define fft_N 256
  4 struct complex
  5 {
  6     double real;
  7     double imag;
  8 };
  9 struct complex c_mul(struct complex a, struct complex b);
 10 void change(struct complex x[]);
 11 struct complex get_WN_kn(double r, double N, double k);
 12 struct complex c_add(struct complex a, struct complex b);
 13 struct complex c_sub(struct complex a, struct complex b);
 14 void fft(struct complex x[], int N);
 15 void ifft(struct complex x[], int N);
 16 void fft1(struct complex x[], struct complex y[]);
 17 double abs2(struct complex a);
 18 double result[fft_N];
 19 struct complex x[fft_N];
 20 unsigned int i, j;
 21 struct complex fft_complex[fft_N];
 22 struct complex fft_result[fft_N];
 23 double fft_double[fft_N];
 24 #include "DSP2833x_Device.h"     // DSP2833x Headerfile Include File
 25 #include "DSP2833x_Examples.h"   // DSP2833x Examples Include File
 26 #include "leds.h"
 27 #include "time.h"
 28 #include "uart.h"
 29 #include "stdio.h"
 30 void I2CA_Init(void);
 31 Uint16 AIC23Write(int Address,int Data);
 32 void Delay(int time);
 33 void Mydelay(Uint32 k);
 34 interrupt void  ISRMcbspSend();
 35 #define LuYin GpioDataRegs.GPADAT.bit.GPIO12
 36 #define LuYin_ST GpioDataRegs.GPADAT.bit.GPIO13
 37 #define BoYin GpioDataRegs.GPADAT.bit.GPIO14
 38 Uint16 temp,i,j,m=0;
 39 Uint32 in_H,in_L;
 40 Uint16 out_H,out_L;
 41 Uint32 in[fft_N];
 42 Uint32 out[fft_N];
 43 void main()
 44 {
 45     InitSysCtrl();
 46     InitPieCtrl();
 47     IER = 0x0000;
 48     IFR = 0x0000;
 49     InitPieVectTable();
 50 
 51     LED_Init();
 52     TIM0_Init(150,200000);//200ms
 53     UARTa_Init(4800);
 54 
 55     InitMcbspaGpio();   //zq
 56     InitI2CGpio();
 57 
 58     I2CA_Init();
 59 
 60     AIC23Write(0x00,0x00);
 61     Delay(100);
 62     AIC23Write(0x02,0x00);
 63     Delay(100);
 64     AIC23Write(0x04,0x7f);
 65     Delay(100);
 66     AIC23Write(0x06,0x7f);
 67     Delay(100);
 68     AIC23Write(0x08,0x14);
 69     Delay(100);
 70     AIC23Write(0x0A,0x00);
 71     Delay(100);
 72     AIC23Write(0x0C,0x00);
 73     Delay(100);
 74     AIC23Write(0x0E,0x43);
 75     Delay(100);
 76     AIC23Write(0x10,0x23);
 77     Delay(100);
 78     AIC23Write(0x12,0x01);
 79     Delay(100);     //AIC23Init
 80 
 81     InitMcbspa();          // Initalize the Mcbsp-A
 82 
 83     EALLOW; // This is needed to write to EALLOW protected registers
 84     PieVectTable.MRINTA = &ISRMcbspSend;
 85     EDIS;   // This is needed to disable write to EALLOW protected registers
 86 
 87     PieCtrlRegs.PIECTRL.bit.ENPIE = 1;   // Enable the PIE block
 88     PieCtrlRegs.PIEIER6.bit.INTx5=1;     // Enable PIE Group 6, INT 5
 89     IER |= M_INT6;                            // Enable CPU INT6
 90 
 91     EINT;   // Enable Global interrupt INTM
 92     ERTM;   // Enable Global realtime interrupt DBGM
 93 
 94     for (i = 0; i <= fft_N - 1; i++)//产生待测序列,模拟在一个周期内采样的点
 95     {
 96         x[i].real = 100* sin(10 * i / (1.0*fft_N)*2.0*PI);
 97         x[i].imag = 0;
 98     }
 99     fft(x, fft_N);
100     for (i = 0; i < fft_N; i++)
101     {
102         result[i] = abs2(x[i]);
103     }
104     while(1)
105     {
106     }
107 }
108 
109 void I2CA_Init(void)
110 {
111    // Initialize I2C
112    I2caRegs.I2CSAR = 0x001A;        // Slave address - EEPROM control code
113 
114    #if (CPU_FRQ_150MHZ)             // Default - For 150MHz SYSCLKOUT
115         I2caRegs.I2CPSC.all = 14;   // Prescaler - need 7-12 Mhz on module clk (150/15 = 10MHz)
116    #endif
117    #if (CPU_FRQ_100MHZ)             // For 100 MHz SYSCLKOUT
118      I2caRegs.I2CPSC.all = 9;       // Prescaler - need 7-12 Mhz on module clk (100/10 = 10MHz)
119    #endif
120 
121    I2caRegs.I2CCLKL = 100;          // NOTE: must be non zero
122    I2caRegs.I2CCLKH = 100;          // NOTE: must be non zero
123    I2caRegs.I2CIER.all = 0x24;      // Enable SCD & ARDY interrupts
124 
125 //   I2caRegs.I2CMDR.all = 0x0020;  // Take I2C out of reset
126    I2caRegs.I2CMDR.all = 0x0420;    // Take I2C out of reset        //zq
127                                     // Stop I2C when suspended
128 
129    I2caRegs.I2CFFTX.all = 0x6000;   // Enable FIFO mode and TXFIFO
130    I2caRegs.I2CFFRX.all = 0x2040;   // Enable RXFIFO, clear RXFFINT,
131 
132    return;
133 }
134 
135 Uint16 AIC23Write(int Address,int Data)
136 {
137 
138 
139    if (I2caRegs.I2CMDR.bit.STP == 1)
140    {
141       return I2C_STP_NOT_READY_ERROR;
142    }
143 
144    // Setup slave address
145    I2caRegs.I2CSAR = 0x1A;
146 
147    // Check if bus busy
148    if (I2caRegs.I2CSTR.bit.BB == 1)
149    {
150       return I2C_BUS_BUSY_ERROR;
151    }
152 
153    // Setup number of bytes to send
154    // MsgBuffer + Address
155    I2caRegs.I2CCNT = 2;
156    I2caRegs.I2CDXR = Address;
157    I2caRegs.I2CDXR = Data;
158    // Send start as master transmitter
159    I2caRegs.I2CMDR.all = 0x6E20;
160    return I2C_SUCCESS;
161 
162 }
163 
164 void Delay(int time)
165 {
166     int i,j,k=0;
167     for(i=0;i<time;i++)
168         for(j=0;j<1024;j++)
169             k++;
170 }
171 
172 void Mydelay(Uint32 k)
173 {
174     while(k--);
175 }
176 
177 interrupt void  ISRMcbspSend(void)
178 {
179     int temp1,temp2;
180 
181     temp1=McbspaRegs.DRR1.all+100;
182     temp2=McbspaRegs.DRR2.all;
183 
184     in_L=McbspaRegs.DRR1.all+100;
185     in_H=McbspaRegs.DRR2.all;
186     for(j=0;j<fft_N-1;j++)
187     {
188         in[j]=in[j+1];
189     }
190     in[fft_N-1]=(in_H<<16)+in_L;
191     for(j=0;j<fft_N;j++)
192     {
193         fft_complex[j].real=in[j];
194         fft_complex[j].imag=0;
195     }
196 
197     if(m >= fft_N)
198     {
199         fft1(fft_complex, fft_result);
200         for(j=0;j<fft_N;j++)
201         {
202             fft_double[j]=abs2(fft_result[j]);
203         }
204         m=0;
205     }
206     m++;
207     McbspaRegs.DXR1.all = temp1;        //放音
208     McbspaRegs.DXR2.all = temp2;
209 
210     PieCtrlRegs.PIEACK.all = 0x0020;
211     //  PieCtrlRegs.PIEIFR6.bit.INTx5 = 0;
212     //  ERTM;
213 }
214 
215 void ifft(struct complex x[], int N)
216 {
217     int i;
218     for (i = 0; i<fft_N; i++)
219     {
220         x[i].imag = -x[i].imag;
221     }
222     fft(x, fft_N);
223 
224     for (i = 0; i<fft_N; i++)
225     {
226         x[i].real = x[i].real*1.0 / fft_N;
227         x[i].imag = -x[i].imag*1.0 / fft_N;
228     }
229 }
230 void fft(struct complex x[], int N)
231 {
232     change(x);
233     struct complex W;
234     struct complex up, down, product;
235     int i, l, j, k;
236     for (i = 0; i<log(fft_N) / log(2); i++)        /*一级蝶形运算 stage */
237     {
238         l = 1 << i;
239         for (j = 0; j<fft_N; j += 2 * l)     /*一组蝶形运算 group,每个group的蝶形因子乘数不同*/
240         {
241             for (k = 0; k<l; k++)        /*一个蝶形运算 每个group内的蝶形运算*/
242             {
243                 W = get_WN_kn(1, fft_N, k * fft_N / (2 * l));
244                 product = c_mul(x[j + k + l], W);
245                 up = c_add(x[j + k], product);
246                 down = c_sub(x[j + k], product);
247                 x[j + k] = up;
248                 x[j + k + l] = down;
249             }
250         }
251     }
252 }
253 void fft1(struct complex x[], struct complex y[])
254 {
255     int n;
256     for (n = 0; n < fft_N; n++)
257     {
258         y[n] = x[n];
259     }
260     change(y);
261     struct complex W;
262     struct complex up, down, product;
263     int i, l, j, k;
264     for (i = 0; i<log(fft_N) / log(2); i++)        /*一级蝶形运算 stage */
265     {
266         l = 1 << i;
267         for (j = 0; j<fft_N; j += 2 * l)     /*一组蝶形运算 group,每个group的蝶形因子乘数不同*/
268         {
269             for (k = 0; k<l; k++)        /*一个蝶形运算 每个group内的蝶形运算*/
270             {
271                 W = get_WN_kn(1, fft_N, k * fft_N / (2 * l));
272                 product = c_mul(y[j + k + l], W);
273                 up = c_add(y[j + k], product);
274                 down = c_sub(y[j + k], product);
275                 y[j + k] = up;
276                 y[j + k + l] = down;
277             }
278         }
279     }
280 }
281 void change(struct complex x[])
282 {
283     int i = 0, j = 0, k = 0;
284     double t;
285     struct complex temp;
286     for (i = 0; i<fft_N; i++)
287     {
288         k = i; j = 0;
289         t = (log(fft_N) / log(2));
290         while ((t--)>0)    //利用按位与以及循环实现码位颠倒
291         {
292             j = j << 1;
293             j |= (k & 1);
294             k = k >> 1;
295         }
296         if (j>i)    //将x(n)的码位互换
297         {
298             temp = x[i];
299             x[i] = x[j];
300             x[j] = temp;
301         }
302     }
303 }
304 struct complex c_mul(struct complex a, struct complex b)//复数乘法
305 {
306     struct complex c;
307     c.real = a.real*b.real - a.imag*b.imag;
308     c.imag = a.real*b.imag + a.imag*b.real;
309     return c;
310 }
311 struct complex c_add(struct complex a, struct complex b)//复数加法
312 {
313     struct complex c;
314     c.imag = a.imag + b.imag;
315     c.real = a.real + b.real;
316     return c;
317 }
318 struct complex c_sub(struct complex a, struct complex b)//复数减法
319 {
320     struct complex c;
321     c.imag = a.imag - b.imag;
322     c.real = a.real - b.real;
323     return c;
324 }
325 struct complex get_WN_kn(double n, double N, double k)//获得旋转因子
326 {
327     int i = 0;
328     struct complex ex;
329     struct complex W;
330     ex.real = 1.0*cos(2 * PI / N);
331     ex.imag = -1.0*sin(2 * PI / N);
332     W.real = 1;
333     W.imag = 0;
334     for (i = 0; i < k*n; i++)
335     {
336         W = c_mul(ex, W);
337     }
338     return W;
339 }
340 double abs2(struct complex a)
341 {
342     double y;
343     y = a.real * a.real + a.imag * a.imag;
344     return y;
345 }

 

标签:频谱,real,struct,音频,complex,3.2,信号,imag
From: https://www.cnblogs.com/Kelvin-Wu/p/16795037.html

相关文章