自相关实验

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

实验要求

  1. 编写一个频率可调、声音大小可调的函数,能产生20Hz到20kHz的响度可调的单频声音信号
  2. 编写一个频率范围可调, 幅度/时长可调的函数, 能产生在20Hz到20KHz的范围内的线性调频信号,信号形式如下:

    $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

  可以听一下原曲:

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

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

  下一个是调频信号,题目中给的公式是一类特定的调频信号,叫作 线性调频(Linear Frequency Modulated)信号,其定义式如下:

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

  其频率会随时间变大而变大,$\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,与实际吻合。

[x,xn]=LFM(1, 0, 500, 0.001, 0, 300); plot(xn,real(x),’b’) [y1, yn1] = sigshift(x, xn, 50); [y2, yn2] = sigshift(x, xn, -100); [y yn] = sigadd(y10.5,yn1,y22,yn2); [z,zn]=myxcorr(y,yn,x,xn); plot(zn,z);