实验一：低通采样定理和内插与抽取实现
一． 实验目的
1. 连续信号和系统的表示方法，以及坊真方法。 2. 用MATLAB 实现连续信号采用与重构的方法， 3. 采样信号的插值和抽取等重采样实现方法。 4. 用时域采样信号重构连续时域信号的原理和方法。
5. 用MATLAB 绘图函数表示信号的基本方法，实验数据的可视化表示。
二． 原理
1 、时域抽样定理
令连续信号xa(t)的傅里叶变换为Xa(jΩ) , 抽样脉冲序列p(t)傅里叶变换为P(jΩ) , 抽样后的信号x^(t)的傅里叶变换为X^(jΩ) 若采用均匀抽样, 抽样周期Ts , 抽样频率为Ωs=2πfs , 由前面分析可知:抽样的过程可以通过抽样脉冲序列p(t)与连续信号xa(t)相乘来完成, 即满足:x^(t)=xa(t) p(t)，又周期信号f(t)傅里叶变换为：
故可以推得p(t)的傅里叶变换为
:
其中
:
根据卷积定理可知
:
得到抽样信号x(t)的傅里叶变换为
:
其表明:信号在时域被抽样后, 他的频谱X(jΩ) 是连续信号频谱X(jΩ) 的形状以抽样频率Ω为间隔周期重复而得到, 在重复过程中幅度被p(t)的傅里叶级数Pn 加权。因为Pn 只是n 的函数, 所以X(jΩ) 在重复的过程中不会使其形状发生变化。
假定信号x(t)的频谱限制在-Ωm~+Ωm 的范围内, 若以间隔Ts 对xa(t)进行抽样, 可知抽样信号X^(t)的频谱X^(jΩ) 是以Ωs 为周期重复。显然, 若在抽样的过程中Ωs
Ωs>=2Ωm 条件,X^(jΩ) 才不会产生频谱的混叠, 接收端完全可以由x^(t)恢复原
连续信号xa(t),这就是低通信号抽样定理的核心内容。 2、信号的重建
从频域看, 设信号最高频率不超过折叠频率:
Xa(jΩ)=Xa(jΩ) |Ω|Ωs/2
则理想取样后的频谱就不会产生混叠, 故有
:
让取样信号x^(t)通过一带宽等于折叠频率的理想低通滤波器:
H(jΩ)=T |Ω|Ωs/2
滤波器只允许通过基带频谱, 即原信号频谱, 故:
Y(jΩ)=X^(jΩ)H(jΩ)=Xa(jΩ)
因此在滤波器的输出得到了恢复的原模拟信号:
y(t)=xa(t)
从时域上看, 上述理想的低通滤波器的脉冲响应为
:
根据卷积公式可求得理想低通滤波器的输出为
:
由上式显然可得
:
则
:
上式表明只要满足取样频率高于两倍信号最高频率, 连续时间函数xa(t)就可用他的取样值xa(nT)来表达而不损失任何信息，这时只要把每一个取样瞬时值与内插函数式相乘求和即可得出xa(t),在每一取样点上, 由于只有该取样值所对应的内插函数式不为零, 所以各个取样点上的信号值不变。
三． 内容
1. 设计连续时间信号线性滤波器分离信号组份
已知信号x (t )=∑(m +1)cos (2π(100m +50) t )，试设计滤波器，分离出如下信号：
m =099
（1）m =1,2,3......50 （2）m =51,52,53......100 （3）m =40,41,42......60 （4）m =1,2,3......40,61,62......100
据以下采样频率：(a) f s =20000Hz (b) f s =10000Hz (c)f s =30000Hz 求信号频谱及相应的滤波器。
参考程序如下：
设计一个Butterworth 模拟带通滤波器，设计指标为：通带频率：1000- 2000Hz ，两侧过渡带宽500Hz ，通带波纹1dB ，阻带衰减100dB 。假设一个信号，其中f1=100Hz,f2=1500Hz,f3=2900Hz。信号的采样频率为10000Hz 。试将原信号与通过该滤波器的模拟信号进行比较。参考程序如下：
wp=[1000 2000]*2*pi;ws=[500 2500]*2*pi;Rp=1;Rs=100; %滤波器设计参数，对于给定Hz 应乘以2
[N,Wn]=buttord(wp,ws,Rp,Rs,'s'); %求得滤波器的最小阶数和截止频率 w=linspace(1,3000,1000)*2*pi; %设置绘制频率响应的频率点 [b,a]=butter(N,Wn,'s'); %设计模拟Butterworth 滤波器
H=freqs(b,a,w); %计算给定频率点的复数频率响应
magH=abs(H);phaH=unwrap(angle(H)); %计算幅频响应和相频响应 plot(w/(2*pi),20*log10(magH)); %以频率为横坐标绘制幅频响应 xlabel('频率/Hz');ylabel('振幅/dB'); title('Butterworth模拟带通滤波器');
hold on;plot([1000 1000],ylim,'r');plot([2000 2000],ylim,'r');%绘带边界grid on figure(2)
dt=1/10000; %模拟信号采样间隔
f1=100;f2=1500;f3=2900;%输入信号的三个频率成分 t=0:dt:0.04; %给定模拟时间段
x=sin(2*pi*f1*t)+0.5*cos(2*pi*f2*t)+0.5*sin(2*pi*f3*t); %输入信号 H=[tf(b,a)]; %滤波器在MATLAB 系统中的表示 [y,t1]=lsim(H,x,t); %模拟输出
subplot(2,1,1),plot(t,x),title('输入信号') %绘出输入信号 subplot(2,1,2),plot(t1,y) %绘制输出信号 title('输出信号'),xlabel('时间/s')
2. 连续时间信号的采样和重建
已知信号x (t )=∑(m +1)cos (2π(100m +50) t )，试以f s =20000Hz 采样频率对该
m =099
信号采样，并用插值公式重建该信号。
参考程序：
①、分别用150HZ 及300HZ 对信号采样
源信号为：
fa=5*sin(2*pi*40*t1)+1.8*sin(4*pi*40*t1)+0.8*sin(5*pi*40*t1)，用150Hz 的频率对f(t)进行采样，其采样图如图1所示；用300Hz 的频率对f(t)进行采样，其采样图如图2所示。程序如下： fs1=150;
t1=-0.1:1/fs1:0.1;
fa=5*sin(2*pi*40*t1)+1.8*sin(4*pi*40*t1)+0.8*sin(5*pi*40*t1); figure(1);plot(t1,fa),xlabel('fs1=150Hz时,fa 采样时域图'); hold off; fs2=300;
t2=-0.1:1/fs2:0.1;
fb=5*sin(2*pi*40*t2)+1.8*sin(4*pi*40*t2)+0.8*sin(5*pi*40*t2); figure(2);plot(t2,fb),xlabel('fs2=300Hz时,fb 采样时域图');
fs1=150Hz时,fa 采样时域图
图1 150HZ 采样频率对信号采样图
fs2=300Hz时,fb 采样时域图
图2 300HZ 采样频率对信号采样图
②、对信号进行快速离散傅里叶变换
将两个采样信号进行快速离散傅里叶变换(FFT),用150Hz 的频率对f(t)进行采样，其采样后快速傅立叶变换频谱图图3所示；用300Hz 的频率对f(t)进行采样，其采样后快速傅立叶变换频谱图图4所示。程序如下：
f=40;fs=150; N=300;k=0:N-1; t=-0.1:1/fs:0.1; w1=150*k/N;
fa=5*sin(2*pi*f*t)+1.8*sin(4*pi*f*t)+0.8*sin(5*pi*f*t); xfa=fft(fa,N);xf1=abs(xfa);
figure(1);plot(w1,xf1),xlabel('fs=150Hz时,fa 经fft 后频谱图. 单位：Hz');
f=40;fs=300; N=300;k=0:N-1; t=-0.1:1/fs:0.1; w2=300*k/N
fb=5*sin(2*pi*f*t)+1.8*sin(4*pi*f*t)+0.8*sin(5*pi*f*t); xfb=fft(fb,N);xf2=abs(xfb);
figure(2);plot(w2,xf2),xlabel('fs=300Hz时,fb 经fft 后频谱图. 单位：Hz ');
fs=150Hz时,fa 经fft 后频谱图. 单位：Hz
图3 150HZ 采样后经FFT 后频谱图
fs=300Hz时,fb 经fft 后频谱图. 单位：Hz
图4 300HZ 采样后经FFT 后频谱图 ③、信号的重建
我们可以通过利用内插法把原信号从采样信号中恢复出来，观察信号在满足怎样的采样条件下能够恢复为原信号，图5和图6分别为恢复后的原信号。程序如下：
Wm=180*pi;Wc=Wm;
fs1=150;Ws=2*pi*fs1;
n=-800:800;nTs1=n/fs1;
fa=5.1*sin(2*pi*40*nTs1)+1.8*sin(4*pi*40*nTs1)+0.8*sin(5*pi*40*nTs1); Dt=1/fs1;t1=-0.1:Dt:0.1;
fa1=fa/fs1*Wc/pi*sinc((Wc/pi)*(ones(length(nTs1),1)*t1-nTs1'*ones(1,length(t1)));
figure(1);plot(t1,fa1); axis([-0.1 0.1 -8 8]);
xlabel('fs=150Hz时,fa 利用内插由样本重建原信号图.'); Wm=180*pi;Wc=Wm; fs2=300;Ws=2*pi*fs2; n=-800:800;nTs2=n/fs2;
fb=5.1*sin(2*pi*40*nTs2)+1.8*sin(4*pi*40*nTs2)+0.8*sin(5*pi*40*nTs2); Dt=1/fs2;t1=-0.1:Dt:0.1;
fb1=fb/fs2*Wc/pi*sinc((Wc/pi)*(ones(length(nTs2),1)*t1-nTs2'*ones(1,length(t1)));
figure(2);plot(t1,fb1); axis([-0.1 0.1 -8 8]);
xlabel('fs=300Hz时,fb 利用内插由样本重建原信号图.'); grid;
fs=150Hz时,fa 利用内插由样本重建原信号图.
图5 150HZ采样后的信号的重建信号
-0.1
-0.08-0.06-0.04-0.0200.020.040.060.080.1
fs=300Hz时,fb 利用内插由样本重建原信号图.
图6 300HZ采样后的信号的重建信号
3. 采样信号的抽取和插值
⑴. 已知信号x (t )=∑(m +1)cos (2π(100m +50) t )，以f s =20000Hz 采样频率采样
m =099
后，设计抽取因子D 和滤波器，分离出如下信号m =1,2,3......50 信号组分。 参考程序如下：
求信号频谱及相应的滤波器。
t = 0:.00025:1; % Time vector x = sin(2*pi*30*t) + sin(2*pi*60*t); figure(1) subplot(211)
stem(t(1:120),x(1:120));hold on y = decimate(x,4);
% View the original and decimated signals: stem(x(1:120)), axis([0 120 -2 2]) % Original signal title('Original Signal') subplot(212)
stem(y(1:30)) % Decimated signal title('Decimated Signal')
⑵. 已知信号x 1(t )=∑(m +1)cos (2π(100m +31) t )，x 2(t )=∑(m +1)cos (2π(100m +37) t ), 以
m =0
m =50
49
99
设计插值因子I 和滤波器，合成信号x 1(t )+x 2(t )。 f s =10000Hz 采样频率采样后，
参考程序如下：
%File_C3:UpSampl.m
%该程序仿真通过零阶保持内插对信号进行上抽样，分析上抽样对信号频谱的影响
clear all clc
f0=0.06; %信号数字频率 N=256;
dt=1; %抽样时间 t=[0:N-1]*dt;
sig=sin(2*pi*f0*t); I=4; %内插因子 N1=N*I;
Addsig(1:N1)=0; for k=1:N
for i=1:I
Addsig(I*k-I+i)=sig(k); %对原始数字信号补零内插，抽样时间变为dL 『I end end
figure(1)
subplot(2,1,1) t2=[0:N1-1]; stem(t,sig); title('Original')
xlabel('Sample time /n') axis([0 19 -1.2 1.2]) subplot(2,1,2)
stem(Addsig);
axis([0 80 -1.2 1.2]) xlabel('Sample time /n')
title('Signal output of Zero-order hold interpolator') Am=fft(sig);
AddAm=fft(Addsig); f1=[0:2/N:2-1/N]; f2=[0:2/N1:2-1/N1]; figure(2)
subplot(2,1,1)
plot(f1,20*log10(abs(Am))); axis([0 1 -20 60]);
ylabel('Amplitude(dB)') xlabel('\omega/ \pi')
title('Amplitude Spectrum') subplot(2,1,2)
plot(f2,20*log10(abs(AddAm))) axis([0 1 -20 60]);
ylabel('Amplitude(dB)') xlabel('\omega/ \pi')
title('Amplitude Spectrum of Zero-order hold interpolation signal') [b,a]=butter(20,1/I); %设计截止频率为pi/l，阶数为10的低通巴特沃思滤波器 y=filter(b,a,Addsig); %对补零后的信号进行低通滤波，完成上抽样过程 UpAm=fft(y); figure(3)
subplot(2,1,1) plot(y);
axis([0 1023 -1.2 1.2])
title('output signal out of Interpolator') xlabel('Sample time /n') subplot(2,1,2)
plot(f2,20*log10(abs(UpAm))) axis([0 1 -30 50]);
ylabel('Amplitude(dB)') xlabel('\omega/ \pi')
%File_C3:UpSampl.m
%该程序仿真通过零阶保持内插对信号进行上抽样，分析上抽样对信号频谱的影响
clear all clc
f0=0.06; %信号数字频率 N=256;
dt=1; %抽样时间
t=[0:N-1]*dt;
sig=sin(2*pi*f0*t);
I=4; %内插因子
N1=N*I;
Addsig(1:N1)=0;
for k=1:N
for i=1:I
Addsig(I*k-I+i)=sig(k); %对原始数字信号补零内插，抽样时间变为dL 『I
end
end
figure(1)
subplot(2,1,1)
t2=[0:N1-1];
stem(t,sig);
title('Original')
xlabel('Sample time /n')
axis([0 19 -1.2 1.2])
subplot(2,1,2)
stem(Addsig);
axis([0 80 -1.2 1.2])
xlabel('Sample time /n')
title('Signal output of Zero-order hold interpolator')
Am=fft(sig);
AddAm=fft(Addsig);
f1=[0:2/N:2-1/N];
f2=[0:2/N1:2-1/N1];
figure(2)
subplot(2,1,1)
plot(f1,20*log10(abs(Am)));
axis([0 1 -20 60]);
ylabel('Amplitude(dB)')
xlabel('\omega/ \pi')
title('Amplitude Spectrum')
subplot(2,1,2)
plot(f2,20*log10(abs(AddAm)))
axis([0 1 -20 60]);
ylabel('Amplitude(dB)')
xlabel('\omega/ \pi')
title('Amplitude Spectrum of Zero-order hold interpolation signal')
[b,a]=butter(20,1/I); %设计截止频率为pi/l，阶数为10的低通巴特沃思滤波器 y=filter(b,a,Addsig); %对补零后的信号进行低通滤波，完成上抽样过程 UpAm=fft(y);
figure(3)
subplot(2,1,1)
plot(y);
axis([0 1023 -1.2 1.2])
title('output signal out of Interpolator') xlabel('Sample time /n')
subplot(2,1,2)
plot(f2,20*log10(abs(UpAm))) axis([0 1 -30 50]);
ylabel('Amplitude(dB)')
xlabel('\omega/ \pi')
四． 设计结果及分析
五． 结论
六． 参考文献
重要参考设计.
理想矩形滤波器的时域表示 clear; clc;
fh=100;
for I=1:4
k=I; fs=k*2*fh;
N=10*k; n=-N:N;
dt=1/fs; T=N*dt;
t=-T:dt:T;
%h=2*fh/fs*sinc(n/k);
%h=sinc(n/k);
subplot(2,2,I)
plot(t,h); hold on；
stem(t,h); hold on;
plot(t,zeros(length(t)),'linewidth',3);
title(['fs/2fh=k,k=',num2str(k)],'fontsize',28); axis('off')
end
实验一：低通采样定理和内插与抽取实现
一． 实验目的
1. 连续信号和系统的表示方法，以及坊真方法。 2. 用MATLAB 实现连续信号采用与重构的方法， 3. 采样信号的插值和抽取等重采样实现方法。 4. 用时域采样信号重构连续时域信号的原理和方法。
5. 用MATLAB 绘图函数表示信号的基本方法，实验数据的可视化表示。
二． 原理
1 、时域抽样定理
令连续信号xa(t)的傅里叶变换为Xa(jΩ) , 抽样脉冲序列p(t)傅里叶变换为P(jΩ) , 抽样后的信号x^(t)的傅里叶变换为X^(jΩ) 若采用均匀抽样, 抽样周期Ts , 抽样频率为Ωs=2πfs , 由前面分析可知:抽样的过程可以通过抽样脉冲序列p(t)与连续信号xa(t)相乘来完成, 即满足:x^(t)=xa(t) p(t)，又周期信号f(t)傅里叶变换为：
故可以推得p(t)的傅里叶变换为
:
其中
:
根据卷积定理可知
:
得到抽样信号x(t)的傅里叶变换为
:
其表明:信号在时域被抽样后, 他的频谱X(jΩ) 是连续信号频谱X(jΩ) 的形状以抽样频率Ω为间隔周期重复而得到, 在重复过程中幅度被p(t)的傅里叶级数Pn 加权。因为Pn 只是n 的函数, 所以X(jΩ) 在重复的过程中不会使其形状发生变化。
假定信号x(t)的频谱限制在-Ωm~+Ωm 的范围内, 若以间隔Ts 对xa(t)进行抽样, 可知抽样信号X^(t)的频谱X^(jΩ) 是以Ωs 为周期重复。显然, 若在抽样的过程中Ωs
Ωs>=2Ωm 条件,X^(jΩ) 才不会产生频谱的混叠, 接收端完全可以由x^(t)恢复原
连续信号xa(t),这就是低通信号抽样定理的核心内容。 2、信号的重建
从频域看, 设信号最高频率不超过折叠频率:
Xa(jΩ)=Xa(jΩ) |Ω|Ωs/2
则理想取样后的频谱就不会产生混叠, 故有
:
让取样信号x^(t)通过一带宽等于折叠频率的理想低通滤波器:
H(jΩ)=T |Ω|Ωs/2
滤波器只允许通过基带频谱, 即原信号频谱, 故:
Y(jΩ)=X^(jΩ)H(jΩ)=Xa(jΩ)
因此在滤波器的输出得到了恢复的原模拟信号:
y(t)=xa(t)
从时域上看, 上述理想的低通滤波器的脉冲响应为
:
根据卷积公式可求得理想低通滤波器的输出为
:
由上式显然可得
:
则
:
上式表明只要满足取样频率高于两倍信号最高频率, 连续时间函数xa(t)就可用他的取样值xa(nT)来表达而不损失任何信息，这时只要把每一个取样瞬时值与内插函数式相乘求和即可得出xa(t),在每一取样点上, 由于只有该取样值所对应的内插函数式不为零, 所以各个取样点上的信号值不变。
三． 内容
1. 设计连续时间信号线性滤波器分离信号组份
已知信号x (t )=∑(m +1)cos (2π(100m +50) t )，试设计滤波器，分离出如下信号：
m =099
（1）m =1,2,3......50 （2）m =51,52,53......100 （3）m =40,41,42......60 （4）m =1,2,3......40,61,62......100
据以下采样频率：(a) f s =20000Hz (b) f s =10000Hz (c)f s =30000Hz 求信号频谱及相应的滤波器。
参考程序如下：
设计一个Butterworth 模拟带通滤波器，设计指标为：通带频率：1000- 2000Hz ，两侧过渡带宽500Hz ，通带波纹1dB ，阻带衰减100dB 。假设一个信号，其中f1=100Hz,f2=1500Hz,f3=2900Hz。信号的采样频率为10000Hz 。试将原信号与通过该滤波器的模拟信号进行比较。参考程序如下：
wp=[1000 2000]*2*pi;ws=[500 2500]*2*pi;Rp=1;Rs=100; %滤波器设计参数，对于给定Hz 应乘以2
[N,Wn]=buttord(wp,ws,Rp,Rs,'s'); %求得滤波器的最小阶数和截止频率 w=linspace(1,3000,1000)*2*pi; %设置绘制频率响应的频率点 [b,a]=butter(N,Wn,'s'); %设计模拟Butterworth 滤波器
H=freqs(b,a,w); %计算给定频率点的复数频率响应
magH=abs(H);phaH=unwrap(angle(H)); %计算幅频响应和相频响应 plot(w/(2*pi),20*log10(magH)); %以频率为横坐标绘制幅频响应 xlabel('频率/Hz');ylabel('振幅/dB'); title('Butterworth模拟带通滤波器');
hold on;plot([1000 1000],ylim,'r');plot([2000 2000],ylim,'r');%绘带边界grid on figure(2)
dt=1/10000; %模拟信号采样间隔
f1=100;f2=1500;f3=2900;%输入信号的三个频率成分 t=0:dt:0.04; %给定模拟时间段
x=sin(2*pi*f1*t)+0.5*cos(2*pi*f2*t)+0.5*sin(2*pi*f3*t); %输入信号 H=[tf(b,a)]; %滤波器在MATLAB 系统中的表示 [y,t1]=lsim(H,x,t); %模拟输出
subplot(2,1,1),plot(t,x),title('输入信号') %绘出输入信号 subplot(2,1,2),plot(t1,y) %绘制输出信号 title('输出信号'),xlabel('时间/s')
2. 连续时间信号的采样和重建
已知信号x (t )=∑(m +1)cos (2π(100m +50) t )，试以f s =20000Hz 采样频率对该
m =099
信号采样，并用插值公式重建该信号。
参考程序：
①、分别用150HZ 及300HZ 对信号采样
源信号为：
fa=5*sin(2*pi*40*t1)+1.8*sin(4*pi*40*t1)+0.8*sin(5*pi*40*t1)，用150Hz 的频率对f(t)进行采样，其采样图如图1所示；用300Hz 的频率对f(t)进行采样，其采样图如图2所示。程序如下： fs1=150;
t1=-0.1:1/fs1:0.1;
fa=5*sin(2*pi*40*t1)+1.8*sin(4*pi*40*t1)+0.8*sin(5*pi*40*t1); figure(1);plot(t1,fa),xlabel('fs1=150Hz时,fa 采样时域图'); hold off; fs2=300;
t2=-0.1:1/fs2:0.1;
fb=5*sin(2*pi*40*t2)+1.8*sin(4*pi*40*t2)+0.8*sin(5*pi*40*t2); figure(2);plot(t2,fb),xlabel('fs2=300Hz时,fb 采样时域图');
fs1=150Hz时,fa 采样时域图
图1 150HZ 采样频率对信号采样图
fs2=300Hz时,fb 采样时域图
图2 300HZ 采样频率对信号采样图
②、对信号进行快速离散傅里叶变换
将两个采样信号进行快速离散傅里叶变换(FFT),用150Hz 的频率对f(t)进行采样，其采样后快速傅立叶变换频谱图图3所示；用300Hz 的频率对f(t)进行采样，其采样后快速傅立叶变换频谱图图4所示。程序如下：
f=40;fs=150; N=300;k=0:N-1; t=-0.1:1/fs:0.1; w1=150*k/N;
fa=5*sin(2*pi*f*t)+1.8*sin(4*pi*f*t)+0.8*sin(5*pi*f*t); xfa=fft(fa,N);xf1=abs(xfa);
figure(1);plot(w1,xf1),xlabel('fs=150Hz时,fa 经fft 后频谱图. 单位：Hz');
f=40;fs=300; N=300;k=0:N-1; t=-0.1:1/fs:0.1; w2=300*k/N
fb=5*sin(2*pi*f*t)+1.8*sin(4*pi*f*t)+0.8*sin(5*pi*f*t); xfb=fft(fb,N);xf2=abs(xfb);
figure(2);plot(w2,xf2),xlabel('fs=300Hz时,fb 经fft 后频谱图. 单位：Hz ');
fs=150Hz时,fa 经fft 后频谱图. 单位：Hz
图3 150HZ 采样后经FFT 后频谱图
fs=300Hz时,fb 经fft 后频谱图. 单位：Hz
图4 300HZ 采样后经FFT 后频谱图 ③、信号的重建
我们可以通过利用内插法把原信号从采样信号中恢复出来，观察信号在满足怎样的采样条件下能够恢复为原信号，图5和图6分别为恢复后的原信号。程序如下：
Wm=180*pi;Wc=Wm;
fs1=150;Ws=2*pi*fs1;
n=-800:800;nTs1=n/fs1;
fa=5.1*sin(2*pi*40*nTs1)+1.8*sin(4*pi*40*nTs1)+0.8*sin(5*pi*40*nTs1); Dt=1/fs1;t1=-0.1:Dt:0.1;
fa1=fa/fs1*Wc/pi*sinc((Wc/pi)*(ones(length(nTs1),1)*t1-nTs1'*ones(1,length(t1)));
figure(1);plot(t1,fa1); axis([-0.1 0.1 -8 8]);
xlabel('fs=150Hz时,fa 利用内插由样本重建原信号图.'); Wm=180*pi;Wc=Wm; fs2=300;Ws=2*pi*fs2; n=-800:800;nTs2=n/fs2;
fb=5.1*sin(2*pi*40*nTs2)+1.8*sin(4*pi*40*nTs2)+0.8*sin(5*pi*40*nTs2); Dt=1/fs2;t1=-0.1:Dt:0.1;
fb1=fb/fs2*Wc/pi*sinc((Wc/pi)*(ones(length(nTs2),1)*t1-nTs2'*ones(1,length(t1)));
figure(2);plot(t1,fb1); axis([-0.1 0.1 -8 8]);
xlabel('fs=300Hz时,fb 利用内插由样本重建原信号图.'); grid;
fs=150Hz时,fa 利用内插由样本重建原信号图.
图5 150HZ采样后的信号的重建信号
-0.1
-0.08-0.06-0.04-0.0200.020.040.060.080.1
fs=300Hz时,fb 利用内插由样本重建原信号图.
图6 300HZ采样后的信号的重建信号
3. 采样信号的抽取和插值
⑴. 已知信号x (t )=∑(m +1)cos (2π(100m +50) t )，以f s =20000Hz 采样频率采样
m =099
后，设计抽取因子D 和滤波器，分离出如下信号m =1,2,3......50 信号组分。 参考程序如下：
求信号频谱及相应的滤波器。
t = 0:.00025:1; % Time vector x = sin(2*pi*30*t) + sin(2*pi*60*t); figure(1) subplot(211)
stem(t(1:120),x(1:120));hold on y = decimate(x,4);
% View the original and decimated signals: stem(x(1:120)), axis([0 120 -2 2]) % Original signal title('Original Signal') subplot(212)
stem(y(1:30)) % Decimated signal title('Decimated Signal')
⑵. 已知信号x 1(t )=∑(m +1)cos (2π(100m +31) t )，x 2(t )=∑(m +1)cos (2π(100m +37) t ), 以
m =0
m =50
49
99
设计插值因子I 和滤波器，合成信号x 1(t )+x 2(t )。 f s =10000Hz 采样频率采样后，
参考程序如下：
%File_C3:UpSampl.m
%该程序仿真通过零阶保持内插对信号进行上抽样，分析上抽样对信号频谱的影响
clear all clc
f0=0.06; %信号数字频率 N=256;
dt=1; %抽样时间 t=[0:N-1]*dt;
sig=sin(2*pi*f0*t); I=4; %内插因子 N1=N*I;
Addsig(1:N1)=0; for k=1:N
for i=1:I
Addsig(I*k-I+i)=sig(k); %对原始数字信号补零内插，抽样时间变为dL 『I end end
figure(1)
subplot(2,1,1) t2=[0:N1-1]; stem(t,sig); title('Original')
xlabel('Sample time /n') axis([0 19 -1.2 1.2]) subplot(2,1,2)
stem(Addsig);
axis([0 80 -1.2 1.2]) xlabel('Sample time /n')
title('Signal output of Zero-order hold interpolator') Am=fft(sig);
AddAm=fft(Addsig); f1=[0:2/N:2-1/N]; f2=[0:2/N1:2-1/N1]; figure(2)
subplot(2,1,1)
plot(f1,20*log10(abs(Am))); axis([0 1 -20 60]);
ylabel('Amplitude(dB)') xlabel('\omega/ \pi')
title('Amplitude Spectrum') subplot(2,1,2)
plot(f2,20*log10(abs(AddAm))) axis([0 1 -20 60]);
ylabel('Amplitude(dB)') xlabel('\omega/ \pi')
title('Amplitude Spectrum of Zero-order hold interpolation signal') [b,a]=butter(20,1/I); %设计截止频率为pi/l，阶数为10的低通巴特沃思滤波器 y=filter(b,a,Addsig); %对补零后的信号进行低通滤波，完成上抽样过程 UpAm=fft(y); figure(3)
subplot(2,1,1) plot(y);
axis([0 1023 -1.2 1.2])
title('output signal out of Interpolator') xlabel('Sample time /n') subplot(2,1,2)
plot(f2,20*log10(abs(UpAm))) axis([0 1 -30 50]);
ylabel('Amplitude(dB)') xlabel('\omega/ \pi')
%File_C3:UpSampl.m
%该程序仿真通过零阶保持内插对信号进行上抽样，分析上抽样对信号频谱的影响
clear all clc
f0=0.06; %信号数字频率 N=256;
dt=1; %抽样时间
t=[0:N-1]*dt;
sig=sin(2*pi*f0*t);
I=4; %内插因子
N1=N*I;
Addsig(1:N1)=0;
for k=1:N
for i=1:I
Addsig(I*k-I+i)=sig(k); %对原始数字信号补零内插，抽样时间变为dL 『I
end
end
figure(1)
subplot(2,1,1)
t2=[0:N1-1];
stem(t,sig);
title('Original')
xlabel('Sample time /n')
axis([0 19 -1.2 1.2])
subplot(2,1,2)
stem(Addsig);
axis([0 80 -1.2 1.2])
xlabel('Sample time /n')
title('Signal output of Zero-order hold interpolator')
Am=fft(sig);
AddAm=fft(Addsig);
f1=[0:2/N:2-1/N];
f2=[0:2/N1:2-1/N1];
figure(2)
subplot(2,1,1)
plot(f1,20*log10(abs(Am)));
axis([0 1 -20 60]);
ylabel('Amplitude(dB)')
xlabel('\omega/ \pi')
title('Amplitude Spectrum')
subplot(2,1,2)
plot(f2,20*log10(abs(AddAm)))
axis([0 1 -20 60]);
ylabel('Amplitude(dB)')
xlabel('\omega/ \pi')
title('Amplitude Spectrum of Zero-order hold interpolation signal')
[b,a]=butter(20,1/I); %设计截止频率为pi/l，阶数为10的低通巴特沃思滤波器 y=filter(b,a,Addsig); %对补零后的信号进行低通滤波，完成上抽样过程 UpAm=fft(y);
figure(3)
subplot(2,1,1)
plot(y);
axis([0 1023 -1.2 1.2])
title('output signal out of Interpolator') xlabel('Sample time /n')
subplot(2,1,2)
plot(f2,20*log10(abs(UpAm))) axis([0 1 -30 50]);
ylabel('Amplitude(dB)')
xlabel('\omega/ \pi')
四． 设计结果及分析
五． 结论
六． 参考文献
重要参考设计.
理想矩形滤波器的时域表示 clear; clc;
fh=100;
for I=1:4
k=I; fs=k*2*fh;
N=10*k; n=-N:N;
dt=1/fs; T=N*dt;
t=-T:dt:T;
%h=2*fh/fs*sinc(n/k);
%h=sinc(n/k);
subplot(2,2,I)
plot(t,h); hold on；
stem(t,h); hold on;
plot(t,zeros(length(t)),'linewidth',3);
title(['fs/2fh=k,k=',num2str(k)],'fontsize',28); axis('off')
end