离散系统的时域分析【数字信号处理二】
离散系统的时域分析
-
一、以下程序中分别使用 conv 和 filter 函数计算 h 和 x 的卷积 y 和 y1，运行程序，并分析 y 和 y1 是否有差别，为什么要使用 x[n]补零后的 x1 来产生 y1；具体分析当 h[n]有 i 个值，x[n]有 j 个值，使用 filter 完成卷积功能，需要如何补零？
-
二、编制程序求解下列两个系统的单位冲激响应和阶跃响应，并绘出其图形。要求分别用 filter、 conv、impz 三种函数完成。给出理论计算结果和程序计算结果并讨论。
-
- 2.1 第一个系统
- 2.2 第二个系统
-
三、听力信号发生器。开发一个具备频率调节功能和音量调节功能的函数模块，在20赫兹至20千赫兹范围内实现可调节音量单频音频信号输出。
-
四、基于自相关算法开发声音信号基频周期检测系统，并分析不同声母及声调条件下的基频特征变化情况。
在所给程序中分别调用conv和filter函数以计算h与x的卷积得到y与y₁。运行该程序后进行结果对比分析：探讨其原因时需考虑x[n]补零后的版本x₁用于生成y₁。详细讨论当h[n]具有i个样本而x[n]具有j个样本时，在filter函数下如何实施卷积运算所需的补零步骤
% Program P2_7
clf;
h = [3 2 1 -2 1 0 -4 0 3]; %impulse response
x = [1 -2 3 -4 3 2 1]; %input sequence
y = conv(h,x);
n = 0:14;
subplot(2,1,1);
stem(n,y);
xlabel('Time index n');
ylabel('Amplitude');
title('Output Obtained by Convolution');
grid;
x1 = [x zeros(1,8)];
y1 = filter(h,1,x1);
subplot(2,1,2);
stem(n,y1);
xlabel('Time index n');
ylabel('Amplitude');
title('Output Generated by Filtering');
grid;
代码解读答：两者之间并无差别；为了实现滤波器处理后的序列与原始序列保持一致，在进行卷积操作时必须确保输入序列经过适当填充才能正确执行卷积运算；具体而言，在完成滤波器作用后得到的结果序列为i + j - 1个样本点；因此在应用滤波器之前需要对原始输入序列x[n]进行前后两端各补充j - 1个零点以满足计算需求。
编制计算机程序以计算下列两个连续时间系统单位冲激响应与阶跃响应，并绘制其冲激响应与阶跃响应的图形。要求分别运用filter、conv、impz三种函数实现。并给出理论推导结果与数值计算结果进行对比分析。
2.1 第一个系统
num=[1,-1];
den=[1,0.75,0.125];
len1=15;
n=0:len1;
n1=0:2*len1;
xn=[1,zeros(1,15)];
yn=ones(1,16);
hn1=filter(num,den,xn);
sn1=filter(num,den,yn);
subplot(2,3,1);
stem(n,hn1);
title("单位冲激响应(filter)");
xlabel("n");ylabel("hn");
subplot(2,3,4);
stem(n,sn1);
title("阶跃响应(filter)");
xlabel("n");ylabel("sn");
hn2=conv(xn,impz(num,den,n));
sn2=conv(yn,impz(num,den,n));
subplot(2,3,2);
stem(n,hn1);
title("单位冲激响应(conv)");
xlabel("n");ylabel("hn");
subplot(2,3,5);
stem(n,sn1);
title("阶跃响应(conv)");
xlabel("n");ylabel("sn");
hn3=impz(num,den,n);
sn3=conv(yn,hn3);
subplot(2,3,3);
stem(n,hn3);
title("单位冲激响应(impz)");
xlabel("n");ylabel("hn");
subplot(2,3,6);
stem(n1,sn3)
title("阶跃响应(impz)");
xlabel("n");ylabel("sn");
axis([0,len1,-0.8, 1]);
代码解读
2.2 第二个系统
num=[0,0.25,0.25,0.25,0.25];
den=1;
len1=30;
n=0:len1;
n1=0:2*len1;
xn=[1,zeros(1,len1)];
yn=ones(1,len1+1);
hn1=filter(num,den,xn);
sn1=filter(num,den,yn);
subplot(2,3,1);
stem(n,hn1);
title("单位冲激响应(filter)");
xlabel("n");ylabel("hn");
subplot(2,3,4);
stem(n,sn1);
title("阶跃响应(filter)");
xlabel("n");ylabel("sn");
hn2=conv(xn,impz(num,den,n));
sn2=conv(yn,impz(num,den,n));
subplot(2,3,2);
stem(n,hn1);
title("单位冲激响应(conv)");
xlabel("n");ylabel("hn");
subplot(2,3,5);
stem(n,sn1);
title("阶跃响应(conv)");
xlabel("n");ylabel("sn");
hn3=impz(num,den,n);
sn3=conv(yn,hn3);
subplot(2,3,3);
stem(n,hn3);
title("单位冲激响应(impz)");
xlabel("n");ylabel("hn");
subplot(2,3,6);
stem(n1,sn3)
title("阶跃响应(impz)");
xlabel("n");ylabel("sn");
axis([0,len1,0,1]);
代码解读
三、听力信号发生器。开发一个具有频率调节功能和音量调节模块的函数模块，并设计并实现一个能够生成20Hz至20kHz范围内的单频声音信号，并确保其响度调节功能的有效性。
function xn=day2ex3()
disp("20~20kHz响度可调单频信号发生器");
T=input("请输入信号时间（秒）：");
fs=44100;
t=0:1/fs:T;
f=input("请输入频率：");
if(f<20||f>20000)
disp("超出频率范围！");
else
a=input("请输入幅度：");
xn=a*sin(2*pi*f*t);
sound(xn,fs)
end
end
代码解读
本节通过自相关算法设计程序来计算声音信号中的基频周期，并对不同声母和声调下的基频特征进行分析。
[voice,fs]=audioread('a_male.wav');
len=length(voice);
%最低周期20ms，取40ms为一帧长度，采样率为16k，N=0.04*16000
N=fs*0.04;
n=10;
dn=0;
temp=voice(N:2*N,1);
subplot(3,1,1);
plot(temp);
title("某一帧");
for i=1:n
temp=voice(N*i+1:N*(i+1),1);
y=xcorr(temp);
[am1,n1]=max(y);
yend=n1-5;
[am2,n2]=max(y(1:yend,1));
dn=dn+n1-n2;
end
d=dn/10;
%时间间隔n/fs ->模拟频率为fs/n
f=fs/d;
subplot(3,1,2);
plot(y);
title("某帧的自相关");
text1 = '采样频率为:%.2fHz\n';
text2 = '基波频率为:%.2fHz\n';
fprintf(text1,fs);
fprintf(text2,f);
f=(-len/2:len/2-1)*fs/len;
X=fftshift(fft(voice));
subplot(3,1,3);
plot(f,abs(X));
title("频谱");
xlabel("f/Hz");
代码解读
采用多帧取均值的方式计算自相关间隔的平均值d=101.4（单位：样本点），其对应的模拟频率估计值则为f=fs/d（其中fs表示采样率）。该基频估计量f₀由自相关法计算得出为157.79 Hz（单位：赫兹），与基于频谱分析所得的结果几乎一致（nearly identical），差异仅为约-0.83 Hz