系列
数字信号处理
2-时域中的离散时间信号
自相关实验

自相关实验

$$ \begin{align*} \newcommand{\dif}{\mathop{}\!\mathrm{d}} \newcommand{\belowarrow}[1]{\mathop{#1}\limits_{\uparrow}} \newcommand{\bd}{\boldsymbol} \newcommand{\L}{\mathscr{L}} \end{align*} $$

讲真我觉得这个实验老师没讲什么有用的。

实验要求

  1. 编写一个频率可调、声音大小可调的函数，能产生20Hz到20kHz的响度可调的单频声音信号

  2. 编写一个频率范围可调, 幅度/时长可调的函数, 能产生在20Hz到20KHz的范围内的线性调频信号，信号形式如下：

    $$ x[nT]=\exp (j2\pi f_c nT + j \pi \mu (nT)^2) $$

    $f_c$ 是中心频率，$T$ 是采样间隔，$\mu$ 是调频率

  3. 将产生的线形调频信号作为输入，并做线性延时相加运算，然后利用互相关对延时信号进行延时估计，观察不同频宽、时长线性调频信号的估计精度差别。

实验过程

  要求1. 很简单，用正弦函数即可，我们要注意的是采样率。matlab 中 sound 函数对采样率的规定如下：

sound(y) % 以默认采样率 2^13=8192 赫兹向扬声器发送音频信号 y。

sound(y,Fs) % 以采样率 Fs 向扬声器发送音频信号 y。

sound(y,Fs,nBits) % 对音频信号 y 使用 nBits 的采样位数。

%附：采样率最小为 1000Hz

  所以容易写出 soundgen（sound generator）：

function y = soundgen(A,f,n,fs)
    y = A*sin(2*pi*f*n/fs);
    sound(y, fs);
end

  我们来产生一个中央C（do~🎵）：

soundgen(1,261.6,1000,1000);

  让我们来唱一下小星星（莫扎特的小星星变奏曲🎵）：

little_star_right=[1 1 5 5 6 6 5 5 4 4 3 3 2 3 1];
little_star_left=[1 8 10 8 11 8 10 8 9 7 8 6 4 5 1];
notes = [1 3 5 6 8 10 12 13 15 17 18 20];
for note = 1:15
    r = 220*2.^((notes(little_star_right(note))+2)./12);
    l = 220*2.^((notes(little_star_left(note))-10)./12);
    soundgen(1,r,1024,2048);
    soundgen(0.3,l,1024,2048);
    pause(0.5)
end

  可以听一下原曲：

{%- include extensions/bilibili.html id='3033530' -%}

  还是差很多，因为钢琴还有一些泛音啥的，所以听起来会好听很多。

  等会儿，不是要弄自相关吗？怎么跑偏了？？？让我们一边听着小星星一边写下面的内容吧！

  下一个是调频信号，题目中给的公式是一类特定的调频信号，叫作 线性调频（Linear Frequency Modulated）信号，其定义式如下：

$$ x(t)=\exp (j2\pi f_c t + j \pi \mu t^2) $$

  通过对相位进行微分，可以求出其频率为：

$$ \frac{1}{2\pi} \frac{\dif \theta}{\dif t} = f_c+\mu t $$

  其频率会随时间变大而变大，$\mu$ 可以看做是“斜率”。

function [x,xn] = LFM(A, fc, u, T, n1, n2)
%generate signal: A*exp(j*2*pi*fc*n*T+j*pi*u*(n*T)^2)
% A  = amplitude
% fc = center frequence
% u  = slope
% T  = sampling time
% n1 = start time
% n2 = end time

    xn=[n1:n2];
    x=A*exp(1i*2*pi*fc*xn*T+1i*pi*u*(xn*T).^2);

end
%为了方便看，相位就从0开始
[x,xn]=LFM(1, 0, 100, 0.001, 0, 300);
plot(xn,real(x),'b')

%为了好玩，可以听一下这个声音
%[x,xn] = LFM(1, 20, 100, 0.001, 0, 3000);
%sound(real(x),1000)

  顺便提一句，由于频率是不断升高的，所以采样频率到后面是不够的，如果后面出现了一些“周期性”的图案，可以调小 $\mu$ 或调大采样频率。

  然后我们对信号进行线性延时相加：（利用之前定义过的 sigadd 和 sigshift，这里再写一次）

function [y,n] = sigadd(x1,n1,x2,n2)
% implements y(n) = x1(n) + x2(n)
% y = sum sequence over n, which includes n1 and n2
% x1 = first sequence over n1
% x2 = second sequence over n2
% 加法器

    n = min(min(n1),min(n2)):max(max(n1),max(n2));
    y1 = zeros(1,length(n));
    y2 = y1;
    y1(find((n>=min(n1))&(n<=max(n1))==1))=x1;
    y2(find((n>=min(n2))&(n<=max(n2))==1))=x2;
    y = y1+y2;

end
function [y,n] = sigshift(x,m,k)
% implement y(n) = x(n-k)
% 延时器

    n = m+k;
    y = x;

end

  我们对信号分别延时 50 和 -100，然后将两个延时信号相加：

[y1, yn1] = sigshift(x, xn, 50);
[y2, yn2] = sigshift(x, xn, -100);
[y yn] = sigadd(y1*0.5,yn1,y2*2,yn2);

  然后终于迎来最重要的自相关啦！matlab 中已经有现成的 xcorr 可用，想要挑战自我的同学也可以根据定义写一个自相关。

function [z,zn] = myxcorr(x,xn,y,yn)
% calculate xcorr(x,y)
    
    %将两个信号补零到相同长度
    n = min(min(xn),min(yn)):max(max(xn),max(yn));
    x1 = zeros(1,length(n));
    y1 = x1;
    x1(find((n>=min(xn))&(n<=max(xn))==1))=x;
    y1(find((n>=min(yn))&(n<=max(yn))==1))=y;
    
    z = xcorr(x1,y1);
    zn = [-length(n)+1:length(n)-1];
end
[z,zn]=myxcorr(y,yn,x,xn);
plot(zn,z);

  从图中可以估算信号分别时移了 -100 和 50，与实际吻合。

时域中的离散时间信号