clear;
f_p=2400; f_s=5000; R_p=3; R_s=25; % 设计要求指标
[n, fn]=buttord(f_p,f_s,R_p,R_s, 's'); % 计算阶数和截止频率
Wn=2*pi*fn; % 转换为角频率
[b,a]=butter(n, Wn, 's'); % 计算H(s)
f=0:100:10000; % 计算频率点和频率范围
s=j*2*pi*f; % s=jw=j*2*pi*f
H_s=polyval(b,s)./polyval(a,s); % 计算相应频率点处H(s)的值
figure(1);
subplot(2,1,1); plot(f, 20*log10(abs(H_s))); % 幅频特性
axis([0 10000 -40 1]);
xlabel('频率 Hz');ylabel('幅度 dB');
subplot(2,1,2); plot(f, angle(H_s)); % 相频特性
xlabel('频率 Hz');ylabel('相角 rad');
figure(2); freqs(b,a); % 也可用指令freqs直接画出H(s)的频率响应曲线。
% ch3example1B.m
clear;
f_p=2400; f_s=5000; R_p=3; R_s=25; % 设计要求指标
[n, fn]=ellipord(f_p,f_s,R_p,R_s,'s'); % 计算阶数和截止频率
Wn=2*pi*fn; % 转换为角频率
[b,a]=ellip(n,R_p,R_s,Wn,'s'); % 计算H(s)
f=0:100:10000; % 计算频率点和频率范围
s=j*2*pi*f; % s=jw=j*2*pi*f
H_s=polyval(b,s)./polyval(a,s); % 计算相应频率点处H(s)的值
figure(1);
subplot(2,1,1); plot(f, 20*log10(abs(H_s))); % 幅频特性
axis([0 10000 -40 1]);
xlabel('频率 Hz');ylabel('幅度 dB');
subplot(2,1,2); plot(f, angle(H_s)); % 相频特性
xlabel('频率 Hz');ylabel('相角 rad');
figure(2); freqs(b,a); % 也可用指令freqs直接画出H(s)的频率响应曲线。
% ch3example2A.m
f_N=8000; % 采样率
f_p=2100; f_s=2500; R_p=3; R_s=25; % 设计要求指标
Ws=f_s/(f_N/2); Wp=f_p/(f_N/2); % 计算归一化频率
[n, Wn]=buttord(Wp,Ws,R_p,R_s); % 计算阶数和截止频率
[b,a]=butter(n, Wn); % 计算H(z)
figure(1);
freqz(b,a, 1000, 8000) % 作出H(z)的幅频相频图, freqz(b,a, 计算点数, 采样率)
subplot(2,1,1); axis([0 4000 -30 3])
figure(2); % 第二种作图方法
f=0:40:4000; % 计算频率点和频率范围
z=exp(j*2*pi*f./(f_N)); %
H_z=polyval(b,z)./polyval(a,z); % 计算相应频率点处H(s)的值
subplot(2,1,1); plot(f, 20*log10(abs(H_z))); % 幅频特性
axis([0 4000 -40 1]);
xlabel('频率 Hz');ylabel('幅度 dB');
subplot(2,1,2); plot(f, angle(H_z)); % 相频特性
xlabel('频率 Hz');ylabel('相角 rad');
% ch3example3A.m
f_N=8000; % 采样率
f_p=1000; f_s=700; R_p=3; R_s=20; % 设计要求指标
Ws=f_s/(f_N/2); Wp=f_p/(f_N/2); % 计算归一化频率
[n, Wn]=cheb1ord(Wp,Ws,R_p,R_s); % 计算阶数和截止频率
[b,a]=cheby1(n, R_p, Wn, 'high'); % 计算H(z)
freqz(b,a, 1000, 8000) % 作出H(z)的幅频相频图, freqz(b,a, 计算点数, 采样率)
subplot(2,1,1); axis([0 4000 -30 3])
% ch3example4A.m
f_N=10000; % 采样率
f_p=[1000, 1500]; f_s=[600, 1900]; R_p=3; R_s=20; % 设计要求指标
Ws=f_s/(f_N/2); Wp=f_p/(f_N/2); % 计算归一化频率
[n, Wn]=ellipord(Wp,Ws,R_p,R_s); % 计算阶数和截止频率
[b,a]=ellip(n, R_p, R_s, Wn); % 计算H(z)
freqz(b,a, 1000, 10000) % 作出H(z)的幅频相频图, freqz(b,a, 计算点数, 采样率)
subplot(2,1,1); axis([0 5000 -30 3])
% ch3example5A.m
f_N=10000; %采样率
f_p=[1000, 1500]; f_s=[1200, 1300]; R_p=3; R_s=30; %设计要求指标
Ws=f_s/(f_N/2); Wp=f_p/(f_N/2); %计算归一化频率
[n, Wn]=cheb2ord(Wp,Ws,R_p,R_s); %计算阶数和截止频率
[b,a]=cheby2(n,R_s, Wn, 'stop'); %计算带阻H(z)系数
freqz(b,a, 1000, 10000) %作出H(z)的幅频相频图, freqz(b,a, 计算点数, 采样率)
subplot(2,1,1); axis([0 5000 -35 3])
% ch3example6A.m
Fs=8000; Ts=1/8000; % 采样率
f0=500; % 梳状滤波器开槽基频率
bw=60/(Fs/2); % 归一化开槽带宽
ab=-3 ; % 计算开槽带宽位置处的衰减分贝值
n=Fs/f0; % 计算滤波器阶数
[num,den] = iircomb(n,bw,ab, 'notch'); % 计算H(z)
freqz(num,den, 4000, 8000); % 作出H(z)的幅频相频图
axis([0 4000 -30 5]);
% ch3example10.m
sim('ch3example10.mdl');
bar([-10:1/5:10-1/5],histdata);
axis([-10 10 0 0.05]);% 作出归一化统计直方图
% ch3example10prg.m
sim('ch3example10.mdl');
bar([-10:1/5:10-1/5],histdata*5);
axis([-10 10 0 0.25]);% 作出归一化统计直方图
% ch3example11prg1.m
t=-1: 0.01 :10 ; % 计算时间范围
f_t=exp(-t).*(t>0);
subplot(3,1,1);plot(t,f_t); % 时域波形
axis([-1 10 -0.1 1.1]); xlabel('time (sec)');
w=-40: 0.1 : 40; % 计算角频率范围
F_w=1./(1+j*w); % 频谱理论结果
subplot(3,1,2);plot(w, abs(F_w)); % 频域幅度谱
axis([-40 40 0 1.1]);xlabel('freq (rad/s)');
subplot(3,1,3);plot(w, angle(F_w)); % 频域相位谱
axis([-40 40 -pi/2 pi/2]);xlabel('freq (rad/s)');
% ch3example11prg2.m
w_m=40; % 截断频率
T=pi/w_m; % 采样间隔
L=5;
t=0: T :L ; % 时域截断
x_t=exp(-t).*(t>0); %信号序列
N=length(x_t); %序列长度(点数)
%-----------------------------------
X_k=fft(x_t); % FFT计算
%-----------------------------------
w0=2*pi/(N*T); % 离散频率间隔
kw=2*pi/(N*T) .*[0: N-1]; % 离散频率样点
X_kw=T.*X_k; % 乘以T得到连续傅利叶变换频谱的样值
%-----------------------------------
subplot(2,1,1);
plot(kw, abs(X_kw),'.','MarkerSize', 14); % 作出数值计算的幅度谱点
axis([-40 90 -0.1 1.1]);xlabel('freq (rad/s)');
hold on; % 保持前面作图曲线不被擦除
w=-40: 0.1 : 40;
X_w=1./(1+j*w); % 理论计算频谱表达式
plot(w, abs(X_w)); % 作图对比
subplot(2,1,2);
plot(kw, angle(X_kw),'.','MarkerSize', 14); % 作出数值计算的相位谱点
hold on;
plot(w, angle(X_w)); % 频域相位谱
axis([-40 90 -pi pi]);xlabel('freq (rad/s)');
% ch3example12prg1.m
Df=5; % 频率间隔
f_s=2*500; % 采样率
N=f_s/Df; % 序列点数
t=0:1./f_s:(N-1)./f_s; % 计算时间段
freq=0:Df:(N-1)*Df; % 计算频率段
f_t=2*sin(2*pi*100*t)+cos(2*pi*180*t); % 信号
F_f=1/f_s*fft(f_t,N); % 用FFT计算频谱
plot(freq-f_s/2,abs(fftshift(F_f))); % 将零频率移动到FFT中心
xlabel('频率 Hz'); ylabel('幅度谱'); % 并作出幅度频谱
% ch3example14prg1.m
fs=200; % 采样率
Delta_f = 1; % 频率分辨率
T = 1/fs; % 时间分辨率
L=1/ Delta_f; % 时域截取长度
N= floor(fs/Delta_f)+1; % 计算截断信号的采样点数
t=0:T:L; % 截取时间段和采样时间点
freq=0: Delta_f :fs; % 分析的频率范围和频率分辨率
f_t =(sin(2*pi*50*t)+0.7*sin(2*pi*75*t))';
% 在截取范围内的分析的信号时域波形
f_t_rectwin= rectwin(N).*f_t./sqrt(sum(abs(rectwin(N).^2))./N);
% 矩形窗（功率归一化）
f_t_hamming= hamming(N) .*f_t./sqrt(sum(abs(hamming(N).^2))./N);
% 海明窗（功率归一化）
f_t_hann = hann(N) .*f_t./sqrt(sum(abs(hann(N).^2))./N);
% 汉宁窗（功率归一化）
F_w_rectwin =T.* fft(f_t_rectwin, N);
% 进行N点FFT，并乘以采样时间间隔T得到频谱
F_w_hamming =T.* fft(f_t_hamming, N); % 加海明窗的频谱
F_w_hann =T.* fft(f_t_hann, N); % 加汉宁窗的频谱figure(1);
subplot(2,2,1);plot(t,f_t);title('Original Signal');
subplot(2,2,2);plot(t, f_t_rectwin);title('Rectwin Windowing');
subplot(2,2,3);plot(t, f_t_hamming);title('hamming Windowing');
subplot(2,2,4);plot(t, f_t_hann);title('hanning Windowing');
figure(2);
subplot(3,1,1);plot(freq, 10*log10(abs(F_w_rectwin)));
title('Rectwin Windowing Spectrum');ylabel('幅度 dB');
axis([0,200,-50,0]);grid on;
subplot(3,1,2);plot(freq, 10*log10(abs(F_w_hamming)));
title('hamming Windowing Spectrum');ylabel('幅度 dB');
axis([0,200,-50,0]);grid on;
subplot(3,1,3);plot(freq, 10*log10(abs(F_w_hann)));
title('hanning Windowing Spectrum');ylabel('幅度 dB');
axis([0,200,-50,0]);grid on;
p_original_signal=var(f_t) % 计算原始信号的平均功率
p_hannwindowed=sum((abs(F_w_hann/T)).^2)/(N^2) % 计算加窗后的信号功率
% ch3example16prg1.m
fs=500; % 采样率
Df=1; % 频率分辨率
N=floor(fs/Df)+1; % 计算的序列点数
t=0:1/fs:(N-1)/fs; % 截取信号的时间段
F=0:Df:fs; % 功率谱估计的频率分辨率和范围
xk=sin(2*pi*50*t)+2*sin(2*pi*130*t)+randn(1,length(t));
% 截取时间段上的离散信号样点序列
Pxx=abs(fft(xk)).^2/(N^2); % 功率谱估计
Pav_timedomaim=sum(xk.^2)/N % 在时域计算信号功率
Pav_freqdomain=sum(Pxx) % 通过功率谱计算信号功率
plot(F,10*log10(Pxx));xlabel('freq Hz');ylabel('PSD dB') % 作出功率谱密度图
% ch3example16prg2.m
fs=500; % 采样率
Df=1; % 频率分辨率
N=floor(fs/Df)+1; % 计算的序列点数
t=0:1/fs:(N-1)/fs; % 截取信号的时间段
F=0:Df:fs; % 功率谱估计的频率分辨率和范围
xk=sin(2*pi*50*t)+2*sin(2*pi*130*t)+randn(1,length(t));
% 截取时间段上的离散信号样点序列
Pxx=(abs(fft(xk(1:167))).^2+...
abs(fft(xk(168:334))).^2+...
abs(fft(xk(335:501))).^2)/3/((N/3)^2);
Pav_timedomaim=sum(xk.^2)/N % 在时域计算信号功率
Pav_freqdomain=sum(Pxx) % 通过功率谱计算信号功率
plot(0:3:fs,10*log10(Pxx));xlabel('freq Hz');ylabel('PSD dB');
%作出功率谱密度图
% ch3example16prg3.m
fs=500; % 采样率
Df=1; % 频率分辨率
N=floor(fs/Df)+1; % 计算的序列点数
t=0:1/fs:(N-1)/fs; % 截取信号的时间段
F=0:Df:fs; % 功率谱估计的频率分辨率和范围
xk=sin(2*pi*50*t)+2*sin(2*pi*130*t)+randn(1,length(t));
% 截取时间段上的离散信号样点序列
Pxx=(abs(fft(xk(1:167))).^2+...
abs(fft(xk(83:249))).^2+...
abs(fft(xk(168:334))).^2+...
abs(fft(xk(250:416))).^2+...
abs(fft(xk(335:501))).^2)/5/((N/3)^2);
Pav_timedomaim=sum(xk.^2)/N % 在时域计算信号功率
Pav_freqdomain=sum(Pxx) % 通过功率谱计算信号功率
plot(0:3:fs,10*log10(Pxx));xlabel('freq Hz');ylabel('PSD dB');
%作出功率谱密度图
% ch3example16prg4.m
fs=500; % 采样率
Df=1; % 频率分辨率
N=floor(fs/Df)+1; % 计算的序列点数
t=0:1/fs:(N-1)/fs; % 截取信号的时间段
F=0:Df:fs; % 功率谱估计的频率分辨率和范围
xk=sin(2*pi*50*t)+2*sin(2*pi*130*t)+randn(1,length(t));
% 截取时间段上的离散信号样点序列
w=hamming(167)';%海明窗
w=w*sqrt(167/sum(w.*w));% 使得海明窗与矩形窗等能量，即加窗后不对信号功率产生影响
Pxx=(abs(fft(w.*xk(1:167))).^2+...
abs(fft(w.*xk(83:249))).^2+...
abs(fft(w.*xk(168:334))).^2+...
abs(fft(w.*xk(250:416))).^2+...
abs(fft(w.*xk(335:501))).^2)/5/((N/3)^2);
Pav_timedomaim=sum(xk.^2)/N % 在时域计算信号功率
Pav_freqdomain=sum(Pxx) % 通过功率谱计算信号功率
plot(0:3:fs,10*log10(Pxx));xlabel('freq Hz');ylabel('PSD dB');%作出功率谱密度图
% ch3example16prg5.m
fs=500; % 采样率
Df=1; % 频率分辨率
N=floor(fs/Df)+1; % 计算的序列点数
t=0:1/fs:(N-1)/fs; % 截取信号的时间段
F=0:Df:fs; % 功率谱估计的频率分辨率和范围
xk=sin(2*pi*50*t)+2*sin(2*pi*130*t)+randn(1,length(t));
% 截取时间段上的离散信号样点序列
[Pxx,F] = psd(xk,512,500,hamming(256),128);
plot(F,10*log10(Pxx/(512/2))); xlabel('freq Hz');ylabel('PSD dB')
Pav=sum(Pxx/(512/2)) % 通过功率谱计算信号功率
% ch3example23prg1.m
kc=150e3; % Hz/V VCO控制灵敏度
omega_n=2*pi*15e3/2.5; % PLL自然角频率
K=2*pi*(0.5*kc); % 估算环路增益
zeta=1; % 临界阻尼
tau_1=K/((omega_n).^2);
tau_2=2*zeta/omega_n-1/K;
freq=0:10:100e3; % 计算频率范围0到100KHz
s=j*2*pi*freq;
G_s=(1+tau_2*s)./(1+tau_1*s); % 环路滤波器传递函数
figure(1);semilogx(freq,(abs(G_s)));% 作出环路滤波器的频率响应
title('环路滤波器幅频响应');xlabel('Hz');ylabel('|G(s)|');
grid on;
b=[tau_2,1]; % 环路滤波器分子系数向量
a=[tau_1,1]; % 环路滤波器分母系数向量
H_s=(G_s*K./s)./(1+G_s*K./s);
figure(2);semilogx(freq,20*log10(abs(H_s)));% 作出闭环频率响应
title('PLL线性相位模型闭环频率响应');xlabel('Hz');ylabel('20log|H(s)|(dB)');
grid on;下载本文
