Chirp信号公式与对离散生成算法之间的差异

讨论产生线性频率变化的公式和它的离散公式之间的差异，并提出Chirp信号的修改方案。

Chirp信号的公式

对于固定频率$f_1$的信号，它的表达式为：$$r\left(t\right)=\cos \left(2\pi \cdot f_1\cdot t\right)$$

其中正弦信号中的相位是频率对时间的积分：$$\theta \left(t\right)=2\pi \cdot {\int }_0^t{f}_1\cdot d\tau =2\pi \cdot f_1\cdot t$$

Chirp信号的频率是变化的，一种最简单的变化就是线性变化，从开始$t=t_0$的频率$f_{start}$，一直线性变化到$t=t_1$时的频率$f_{end}$。那么任何时刻的频率值为：$$f\left(t\right)=\frac{f_{end}-f_{start}}{t_1-t_0}\cdot t+f_{start}$$

对应的正弦信号相角为：$$\theta \left(t\right)=2\pi \cdot {\int }_{t_0}^{t_1}f\left(\tau \right)d\tau =2\pi {\int }_{t_0}^{t_1}\frac{f_{end}-f_{start}}{t_1-t_0}\cdot \tau +f_{start}d\tau$$$$=2\pi \left[\frac{f_{end}-f_{start}}{t_1-t_0}\cdot \frac{1}{2}\left(t_{}^2-t_0^2\right)+f_{start}\left(t_{}-t_0\right)\right]$$

如果假设$t_0=0$，那么上面的公式就是：$$\theta \left(t\right)=2\pi \left[\frac{f_{end}-f_{start}}{t_1}\cdot \frac{1}{2}{t}^2+f_{start}\cdot t\right]$$

所以，Chirp信号的公式为：
$$r\left(t\right)=\cos \left[2\pi \cdot \left(\frac{f_{end}-f_{start}}{t_1}\cdot \frac{1}{2}{t}^2+f_{start}\cdot t\right)\right]$$

下面是在智能车比赛中相应的参数：

起始时间：$t_0=0s$
结束时间：$t_1=0.2048s$
起始频率： $f_{start}=250Hz$
结束频率：$f_{end}=2000Hz$

利用时间间隔$t_s=1/f_s=100\mu s$对该信号进行离散化，得到的波形为：

def chirpf(t):
    return cos(2*pi*((fend-fstart)/t1*0.5*t*t+fstart*t))

  

 

  
1
  
2
 

使用MATLAB产生的Chirp信号波形：

tend=0.2048
fs = 1e4
fstart=250
fend=2000
tdim = 0:1/fs:(tend-1/fs)
chirpdim=chirp(tdim,fstart,tend,fend)

  

 

  
1
  
2
  
3
  
4
  
5
  
6
 

import scipy.signal
data = scipy.signal.chirp(tdim, fstart,t1, fend)

  

 

  
1
  
2
 

通过数值对比，可以看到上面三种方式得到的Chirp信号的误差为0。下图显示了公式所计算得到的Chirp信号与MATLAB所产生的Chirp信号的差值。其中微小的误差可能来源于MATLAB结果经过clipboard剪切板之后，数值显示精度降低所引起的。

另外一部分则可能是对于数据的的最后一点的频率微小差异的定义。

Chirp信号离散生成算法

下面是智能车竞赛信号板单片机通过DAC输出Chirp信号对输出信号buffer进行初始化的程序。

//------------------------------------------------------------------------------
void InitDACBuffer(float fStartF, float fEndF) {
    float fAngle = 0;
    float fFrequency;
    float fDeltaT = 1.0 / DAC_OUTPUT_FREQUENCY;
    
    int i;
    for(i = 0; i < DAC_BUFFER; i ++) {        
        g_nDACBuffer[i] = (unsigned short)((sin(fAngle * 2 * 3.1415926) + 1.0) / 2 * 0x4ff)+0x100;
        fFrequency = (fEndF - fStartF) * (i + 1) / DAC_BUFFER + fStartF;
        fAngle += fFrequency * fDeltaT;
    }
}

  

 

  
1
  
2
  
3
  
4
  
5
  
6
  
7
  
8
  
9
  
10
  
11
  
12
  
13
 

在程序初始化期间调用如下命令：

    InitDACBuffer(250, 2000);

  

 

  
1
 

根据上面产生的程序，下面绘制出对应的波形：

DAC_BUFFER = 2048
dacbuf = []
angle = 0
deltat = 1/1e4

for i in range(0, DAC_BUFFER):
    data = cos(angle*2*pi)
    dacbuf.append(data)
    frequency = (fend - fstart) * (i + 1) / DAC_BUFFER + fstart
    angle = angle + frequency * deltat

plt.plot(tdim, dacbuf)
plt.xlabel('Time(s)')
plt.ylabel('Amplitude')
plt.grid(True)
plt.show()

  

 

  
1
  
2
  
3
  
4
  
5
  
6
  
7
  
8
  
9
  
10
  
11
  
12
  
13
  
14
  
15
  
16
 

对比由离散程序所产生的数据与前面通过公式计算出的波形，可以看到它们之间存在着差异。下面将公式所计算出的数值与前面程序所产生信号相加，所得到的误差信号：

使用程序产生Chirp信号的时候，信号的相位是通过离散积分而得到的，其中所出现的积分误差随着时间增加而叠加，最终使得所产生的信号的相位与实际的相位越来越大，信号的差别，也就增加了。

修改方案

将最终信号Chirp生成的公式修改为：

$$t_n=n\cdot t_s=\frac{n}{f_s}=n\cdot 10^{-4}$$
$$x\left[n\right]=\sin \left[2\pi \cdot \left(\frac{f_{end}-f_{start}}{t_1}\cdot \frac{1}{2}{t}_n^2+f_{start}\cdot t_n\right)\right]$$

使用sin信号的目的，就是防止最初的信号开始的跳变。根据上面公式所产生的波形为：

生成12bit DAC转换对应的数据：

for i in range(DAC_BUFFER):
    tn = i / fs
    angle = (fend-fstart)/t1 * tn * tn / 2 + fstart * tn
    datasin = sin(2*pi*angle)

    dataint16 = int((datasin + 1.0) / 2 * 0x4ff + 0x100)
    dacbuf.append(dataint16)

  

 

  
1
  
2
  
3
  
4
  
5
  
6
  
7
 

信号波形为：

生成相关的整型信号，便于MCU输出模拟信号：

#------------------------------------------------------------
fstart = 250
fend = 2000
t0 = 0
t1 = 0.2048
ts = 100e-6

tdim = linspace(t0, t1, int((t1-t0)/ts), endpoint=False)

#------------------------------------------------------------
DAC_BUFFER = 2048
dacbuf = []
fs = 1e4

for i in range(DAC_BUFFER):
    tn = i / fs
    angle = (fend-fstart)/t1 * tn * tn / 2 + fstart * tn
    datasin = sin(2*pi*angle)

    dataint16 = int((datasin + 1.0) / 2 * 0x3f)
    dacbuf.append(dataint16)

tspsave('chirp6bit', dacbuf=dacbuf)

  

 

  
1
  
2
  
3
  
4
  
5
  
6
  
7
  
8
  
9
  
10
  
11
  
12
  
13
  
14
  
15
  
16
  
17
  
18
  
19
  
20
  
21
  
22
  
23
 

下面是将NPZ数据转换成C51的变量的程序。

dacbuf0 = tspload('chirp8bit', 'dacbuf')
dacbuf1 = tspload('chirp7bit', 'dacbuf')
dacbuf2 = tspload('chirp6bit', 'dacbuf')

pastestr = 'unsigned char const code g_ucChirpData[] = {\r\n'
linenumber = 16
lines = int(len(dacbuf2) / linenumber)

for i in range(lines):
    linedata = 127 - dacbuf2[i * linenumber : (i + 1) * linenumber]

    datastr = str(linedata).strip('[]').split(' ')
    datastr1 = [s for s in datastr if len(s) > 0]
    linestr = ','.join(datastr1).lstrip(',') + ',\r\n'



    pastestr = pastestr + '    ' + linestr

pastestr = pastestr + '};\r\n'

clipboard.copy(pastestr)
printf(pastestr)

  

 

  
1
  
2
  
3
  
4
  
5
  
6
  
7
  
8
  
9
  
10
  
11
  
12
  
13
  
14
  
15
  
16
  
17
  
18
  
19
  
20
  
21
  
22
  
23
 

文章来源: zhuoqing.blog.csdn.net，作者：卓晴，版权归原作者所有，如需转载，请联系作者。

原文链接：zhuoqing.blog.csdn.net/article/details/105762739