11 实验目的
1．了解数字信号处理系统的一般构成；
2．掌握奈奎斯特抽样定理。
12 实验仪器
1．YBLD智能综合信号源测试仪 1台
2．双踪示波器 1台
3．MCOM－TG305数字信号处理与现代通信技术实验箱 1台
4．PC机（装有MATLAB、MCOM－TG305配套实验软件） 1台
13 实验原理
一个典型的DSP系统除了数字信号处理部分外，还包括A/D和D/A两部分。这是因为自然界的信号，如声音、图像等大多是模拟信号，因此需要将其数字化后进行数字信号处理，模拟信号的数字化即称为A/D转换。数字信号处理后的数据可能需还原为模拟信号，这就需要进行D/A转换。一个仅包括A/D和D/A两部分的简化数字信号处理系统功能如图1所示。
A/D转换包括三个紧密相关的过程，即抽样、量化和编码。A/D转换中需解决的以下几个重要问题：抽样后输出信号中还有没有原始信号的信息？如果有能不能把它取出来？抽样频率应该如何选择？
奈奎斯特抽样定理（即低通信号的均匀抽样定理）告诉我们，一个频带限制在0至fx以内的低通信号x(t)，如果以fs≥2fx的抽样速率进行均匀抽样，则x(t)可以由抽样后的信号xs(t)完全地确定，即xs(t)包含有x(t)的成分，可以通过适当的低通滤波器不失真地恢复出x(t)。最小抽样速率fs=2fx称为奈奎斯特速率。
低通
译码
编码
量化
抽样
输入信号 样点输出 滤波输出
A/D（模数转换） D/A（数模转换）
图1 低通采样定理演示
为方便实现，实验中更换了一种表现形式，即抽样频率固定（10KHz），通过改变输入模拟信号的频率来展示低通抽样定理。我们可以通过研究抽样频率和模拟信号最高频率分量的频率之间的关系，来验证低通抽样定理。
14 实验内容
1．软件仿真实验：编写并调试MATLAB程序，分析有关参数，记录有关波形。
2．硬件实验：输入不同频率的正弦信号，观察采样时钟波形、输入信号波形、样点输出波形和滤波输出波形。
15 MATLAB参考程序和仿真内容
%%
%f—余弦信号的频率
% M—基2 FFT幂次数 N=2^M为采样点数，这样取值是为了便于作基2的FFT分析
%2 采样频率Fs
%%
function samples(f,Fs,M)
N=2^M; % fft点数=取样总点数
Ts=1/Fs; % 取样时间间隔
T=NTs; % 取样总时间=取样总点数取样时间间隔
n=0:N-1;
t=nTs;
Xn=cos(2fpit);
subplot(2,1,1);
stem(t,Xn);
axis([0 T 11min(Xn) 11max(Xn)]);
xlabel('t -->');
ylabel('Xn');
Xk=abs(fft(Xn,N));
subplot(2,1,2);
stem(n,Xk);
axis([0 N 11min(Xk) 11max(Xk)]);
xlabel('frequency -->');
ylabel('!Xk!');
%%
假如有一个1Hz的余弦信号y=cos(2πt)，对其用4Hz的采样频率进行采样，共采样32点，只需执行samples(1,4,5)，即可得到仿真结果。
软件仿真实验内容如下表所示：
仿真参数
f
Fs
Wo(计算)
Xn（图形）
Xk（图形）
（1，4，5）
另外记录图形，并标图号
（1，8，5）
（2，8，6）
自 选
16 硬件实验步骤
本实验箱采样频率fs固定为10KHz，低通滤波器的截止频率约为45KHz。
1、用低频信号源产生正弦信号，正弦信号源频率f自定，并将其接至2TP2（模拟输入）端，将示波器通道一探头接至2TP6（采样时钟）端观察采样时钟波形，示波器通道二探头接至2TP2观察并记录输入信号波形。
2、将示波器通道二探头接至2TP3观察并记录样点输出波形。
3、将示波器通道二探头接至2TP4观察并记录滤波输出波形。
4、根据采样定理，分f=fs /8、f=fs/4、f=fs/2等3种情况更改正弦信号频率，重复步骤2至步骤3。
5、用低频信号源产生方波信号，重复步骤1至步骤4。
17 思考题
1、 讨论在仿真实验中所计算的数字域频率Wo和Xk的图形中非零谱线位置之间的对应关系。
2、 讨论在仿真实验中自选参数的意义。
3、将在2TP2端加方波信号后的恢复波形，与相同频率的正弦信号的恢复波形相比，能够得出哪些结论？
2 FFT频谱分析实验
21 实验目的
1．通过实验加深对快速傅立叶变换（FFT）基本原理的理解。
2．了解FFT点数与频谱分辨率的关系，以及两种加长序列FFT与原序列FFT的关系。
22 实验仪器
1．YBLD智能综合信号源测试仪 1台
2．双踪示波器 1台
3．MCOM－TG305数字信号处理与现代通信技术实验箱 1台
4．PC机（装有MATLAB、MCOM－TG305配套实验软件） 1台
23 实验原理
离散傅里叶变换（DFT）和卷积是信号处理中两个最基本也是最常用的运算，它们涉及到信号与系统的分析与综合这一广泛的信号处理领域。实际上卷积与DFT之间有着互通的联系：卷积可化为DFT来实现，其它的许多算法，如相关、滤波和谱估计等都可化为DFT来实现，DFT也可化为卷积来实现。
对N点序列x(n)，其DFT变换对定义为：
在DFT运算中包含大量的重复运算。FFT算法利用了蝶形因子WN的周期性和对称性，从而加快了运算的速度。FFT算法将长序列的DFT分解为短序列的DFT。N点的DFT先分解为2个N/2点的DFT，每个N/2点的DFT又分解为2个N/4点的DFT。按照此规律，最小变换的点数即所谓的“基数（radix）。”因此，基数为2的FFT算法的最小变换（或称蝶形）是2点DFT。一般地，对N点FFT，对应于N个输入样值，有N个频域样值与之对应。一般而言，FFT算法可以分为时间抽取（DIT）FFT和频率抽取（DIF）两大类。
在实际计算中，可以采用在原来序列后面补0的加长方法来提高FFT的分辨率；可以采用在原来序列后面重复的加长方法来增加FFT的幅度。
24 实验内容
1．软件仿真实验：分别观察并记录正弦序列、方波序列及改变FFT的点数后的频谱；分别观察并记录正弦序列、方波序列及2种加长序列等信号的频谱。
2．硬件实验：分别观察并记录正弦信号、方波信号及改变FFT的点数后的频谱。
25 MATLAB参考程序和仿真内容
%%
function[x]=ffts(mode,M)
Nfft=2^M;
x=zeros(1,Nfft); %定义一个长度为Nfft的一维全0数组
if mode= =1 for n=0:Nfft-1 x(n+1)=sin(2pin/Nfft); end
end %定义一个长度为Nfft的单周期正弦序列
if mode= =2 for n=0:Nfft-1 x(n+1)=sin(4pin/Nfft); end
end %定义一个长度为Nfft的双周期正弦序列
if mode= =3 for n=0:Nfft/2-1 x(n+1)=sin(4pin/Nfft); end
end %定义一个长度为Nfft/2的正弦序列，后面一半为0序列。
if mode= =4 for n=0:Nfft-1 x(n+1)=square(2pin/Nfft); end
end
if mode= =5 for n=0:Nfft-1 x(n+1)=square(2pin/Nfft); end
end
if mode= =6 for n=0:Nfft/2-1 x(n+1)=square(4pin/Nfft); end
end
n=0:Nfft-1;
subplot(2,1,1);
stem(n,x);
axis([0 Nfft-1 11min(x) 11max(x)]);
xlabel('Points-->');
ylabel('x(n)');
X=abs(fft(x,Nfft));
subplot(2,1,2);
stem(n,X);
axis([0 Nfft-1 11min(X) 11max(X)]);
xlabel('frequency-->');
ylabel('!X(k)!');
%%
假设需观察方波信号的频谱，对一个周期的方波信号作32点的FFT，则只需在MATLAB的命令窗口下键入：[x]=ffts(21,5) ，程序进行模拟，并且输出FFT的结果。
关于软件仿真实验内容，建议在完成大量仿真例子的基础上，选择能够体现实验要求的4个以上的例子进行记录。例如要观察后面补0的加长方法来提高FFT的分辨率的现象，可以仿真ffts(4,5)和ffts(6,6)两个例子。
26 硬件实验步骤
1．将低频信号源输出加到实验箱模拟通道1输入端，将示波器探头接至模拟通道1输出端。
2．在保证实验箱正确加电且串口电缆连接正常的情况下，运行数字信号处理与DSP应用实验开发软件，在“数字信号处理实验”菜单下选择“FFT频谱分析”子菜单，出现显示FFT频谱分析功能提示信息的窗口。
3．用低频信号产生器产生一个1KHz的正弦信号。
4．选择FFT频谱分析与显示的点数为64点，开始进行FFT运算。此后，计算机将周期性地取回DSP运算后的FFT数据并绘图显示
5．改信号源频率，观察并记录频谱图的变化。
6．选择FFT的点数为128点，观察并记录频谱图的变化。
7．更改正弦信号的频率，重复步骤4 ~步骤6。
8．用低频信号产生器产生一个1KHz的方波信号，重复步骤4 ~步骤7。注意：应根据实验箱采样频率fs为10KHz和方波信号的频带宽度选择方波信号的频率。
本硬件实验要进行两种信号，每个信号两种频率，每个信号两种点数等共8次具体实验内容，性质能够体现实验要求的4个以上的例子进行记录。
27 思考题
1．对同一个信号，不同点数FFT观察到的频谱图有何区别？
2．序列加长后FFT与原序列FFT的关系是什么，试推导其中一种关系。
3．用傅立叶级数理论，试说明正弦信号频谱和方波信号频谱之间的关系。
3 IIR滤波器设计实验
31 实验目的
1．通过实验加深对IIR滤波器基本原理的理解。
2．学习编写IIR滤波器的MATLAB仿真程序。
32 实验仪器
1．YBLD智能综合信号源测试仪 1台
2．双踪示波器 1台
3．MCOM－TG305数字信号处理与现代通信技术实验箱 1台
4．PC机（装有MATLAB、MCOM－TG305配套实验软件） 1台
33 实验原理
IIR滤波器有以下几个特点：
1．IIR数字滤波器的系统函数可以写成封闭函数的形式。
2．IIR数字滤波器采用递归型结构，即结构上带有反馈环路。IIR滤波器运算结构通常由延时、乘以系数和相加等基本运算组成，可以组合成直接型、正准型、级联型、并联型四种结构形式，都具有反馈回路。由于运算中的舍入处理，使误差不断累积，有时会产生微弱的寄生振荡。
3．IIR数字滤波器在设计上可以借助成熟的模拟滤波器的成果，如巴特沃斯、契比雪夫和椭圆滤波器等，有现成的设计数据或图表可查，其设计工作量比较小，对计算工具的要求不高。在设计一个IIR数字滤波器时，我们根据指标先写出模拟滤波器的公式，然后通过一定的变换，将模拟滤波器的公式转换成数字滤波器的公式。
4．IIR数字滤波器的相位特性不好控制，对相位要求较高时，需加相位校准网络。
在MATLAB下设计IIR滤波器可使用Butterworth函数设计出巴特沃斯滤波器，使用Cheby1函数设计出契比雪夫I型滤波器，使用Cheby2设计出契比雪夫II型滤波器，使用ellipord函数设计出椭圆滤波器。下面主要介绍前两个函数的使用。
与FIR滤波器的设计不同，IIR滤波器设计时的阶数不是由设计者指定，而是根据设计者输入的各个滤波器参数（截止频率、通带滤纹、阻带衰减等），由软件设计出满足这些参数的最低滤波器阶数。在MATLAB下设计不同类型IIR滤波器均有与之对应的函数用于阶数的选择。
一、巴特沃斯IIR滤波器的设计
在MATLAB下，设计巴特沃斯IIR滤波器可使用butter函数。
Butter函数可设计低通、高通、带通和带阻的数字和模拟IIR滤波器，其特性为使通带内的幅度响应最大限度地平坦，但同时损失截止频率处的下降斜度。在期望通带平滑的情况下，可使用butter函数。
butter函数的用法为：
[b,a]=butter(n,Wn,/ftype/)
其中n代表滤波器阶数，Wn代表滤波器的截止频率，这两个参数可使用buttord函数来确定。buttord函数可在给定滤波器性能的情况下，求出巴特沃斯滤波器的最小阶数n，同时给出对应的截止频率Wn。buttord函数的用法为：
[n,Wn]= buttord(Wp,Ws,Rp,Rs)
其中Wp和Ws分别是通带和阻带的拐角频率（截止频率），其取值范围为0至1之间。当其值为1时代表采样频率的一半。Rp和Rs分别是通带和阻带区的波纹系数。
不同类型（高通、低通、带通和带阻）滤波器对应的Wp和Ws值遵循以下规则：
1．高通滤波器：Wp和Ws为一元矢量且Wp>Ws；
2．低通滤波器：Wp和Ws为一元矢量且Wp<Ws；
3．带通滤波器：Wp和Ws为二元矢量且Wp<Ws，如Wp=[02,07],Ws=[01,08];
4．带阻滤波器：Wp和Ws为二元矢量且Wp>Ws，如Wp=[01,08],Ws=[02,07]。
二、契比雪夫I型IIR滤波器的设计
在期望通带下降斜率大的场合，应使用椭圆滤波器或契比雪夫滤波器。在MATLAB下可使用cheby1函数设计出契比雪夫I型IIR滤波器。
cheby1函数可设计低通、高通、带通和带阻契比雪夫I型滤IIR波器，其通带内为等波纹，阻带内为单调。契比雪夫I型的下降斜度比II型大，但其代价是通带内波纹较大。
cheby1函数的用法为：
[b,a]=cheby1(n,Rp,Wn,/ftype/)
在使用cheby1函数设计IIR滤波器之前，可使用cheblord函数求出滤波器阶数n和截止频率Wn。cheblord函数可在给定滤波器性能的情况下，选择契比雪夫I型滤波器的最小阶和截止频率Wn。
cheblord函数的用法为：
[n,Wn]=cheblord(Wp,Ws,Rp,Rs)
其中Wp和Ws分别是通带和阻带的拐角频率（截止频率），其取值范围为0至1之间。当其值为1时代表采样频率的一半。Rp和Rs分别是通带和阻带区的波纹系数。
34 实验内容
1．软件仿真实验：编写并调试MATLAB程序，选择不同形式，不同类型的4种滤波器进行仿真，记录幅频和相频特性，对比巴特沃斯滤波器和契比雪夫滤波器。
2．硬件实验：设计IIR滤波器，在计算机上观察冲激响应、幅频特性和相频特性，然后下载到实验箱。用示波器观察输入输出波形，测试滤波器的幅频响应特性。
35 MATLAB参考程序和仿真内容
%%
%mode: 1--巴特沃斯低通；2--巴特沃斯高通；3--巴特沃斯带通；4--巴特沃斯带阻
% 5--契比雪夫低通；6--契比雪夫高通；7--契比雪夫带通；8--契比雪夫带阻
%fp1,fp2： 通带截止频率，当高通或低通时只有fp1有效
%fs1, fs2： 阻带截止频率，当高通或低通时只有fs1有效
%rp: 通带波纹系数
%as: 阻带衰减系数
%sample: 采样率
%h: 返回设计好的滤波器系数
%%
function[b,a]=iirfilt(mode,fp1,fp2,fs1,fs2,rp,as,sample)
wp1=2fp1/sample;wp2=2fp2/sample;
ws1=2fs1/sample;ws2=2fs2/sample;
%得到巴特沃斯滤波器的最小阶数N和3bd频率wn
if mode<3[N,wn]=buttord(wp1,ws1,rp,as);
elseif mode<5[N,wn]=buttord([wp1 wp2],[ws1 ws2],rp,as);
%得到契比雪夫滤波器的最小阶数N和3bd频率wn
elseif mode<7[N,wn]=cheb1ord(wp1,ws1,rp,as);
else[N,wn]=cheblord([wp1 wp2],[ws1 ws2],rp,as);
end
%得到滤波器系数的分子b和分母a
if mode= =1[b,a]=butter(N,wn);end
if mode= =2[b,a]=butter(N,wn,/high/);end
if mode= =3[b,a]=butter(N,wn);end
if mode= =4[b,a]=butter(N,wn,/stop/);end
if mode= =5[b,a]=cheby1(N,rp,wn);end
if mode= =6[b,a]=cheby1(N,rp,wn,/high/);end
if mode= =7[b,a]=cheby1(N,rp,wn);end
if mode= =8[b,a]=cheby1(N,rp,wn,/stop/);end
set(gcf,/menubar/,menubar);
freq_response=freqz(b,a);
magnitude=20log10(abs(freq_response));
m=0:511;
f=msample/(2511);
subplot(3,1,1);plot(f,magnitude);grid; %幅频特性
axis([0 sample/2 11min(magnitude) 11max(magnitude)]);
ylabel('Magnitude');xlabel('Frequency-->');
phase=angle(freq_response);
subplot(3,1,2);plot(f,phase);grid; %相频特性
axis([0 sample/2 11min(phase) 11max(phase)]);
ylabel('Phase');xlabel('Frequency-->');
h=impz(b,a,32); %32点的单位函数响应
t=1:32;
subplot(3,1,3);stem(t,h);grid;
axis([0 32 12min(h) 11max(h)]);
ylabel('h(n)');xlabel('n-->');
%%
假设需设计一个巴特沃斯低通IIR滤波器，通带截止频率为2KHz，阻带截止频率为3KHz，通带波纹系数为1，阻带衰减系数为20，采样频率为10KHz，则只需在MATLAB的命令窗口下键入：
[b,a]=iirfilt(1,2000,3000,2400,2600,1,20,10000)
程序进行模拟，并且按照如下顺序输出数字滤波器系统函数
的系数
b= b0 b1 ……bn
a= a0 a1 ……an
关于软件仿真实验内容，建议在完成大量仿真例子的基础上，选择能够体现实验要求的4个例子进行记录，系统函数只要记录系统的阶数。
36 硬件实验步骤
1．根据实验箱采样频率fs为10KHz的条件，用低频信号发生器产生一个频率合适的低频正弦信号，将其加到实验箱模拟通道1输入端，将示波器通道1探头接至模拟通道1输入端，通道2探头接至模拟通道2输出端。
2．在保证实验箱正确加电且串口电缆连接正常的情况下，运行数字信号处理与DSP应用实验开发软件，在“数字信号处理实验”菜单下选择“IIR滤波器”子菜单，出现提示信息。
3．输入滤波器类型、滤波器截止频率等参数后，分别点击“幅频特性”和“相频特性”按钮，在窗口右侧观察IIR滤波器的幅频特性和相频特性。此时提示信息将消失，如需查看提示信息，可点击“设计说明”按钮。
4．点击“下载实现”按钮，IIR滤波器开始工作，此时窗口右侧将显示IIR滤波器的幅频特性。
5．根据输入滤波器类型，更改低频信号源的频率，观察示波器上输入输出波形幅度的变化情况，测量IIR滤波器的幅频响应特性，看其是否与设计的幅频特性一致。
6．更改滤波器类型、滤波器截止频率等参数（共4种），重复步骤3至步骤5。所选择的例子参数最好和MATLAB仿真程序的例子一样。
7．用低频信号产生器产生一个500Hz的方波信号，分别设计3种滤波器，完成如下表要求的功能，并且记录参数和波形。
功 能
滤波器类型
参 数
输出波形
fp1
fp2
fs1
fs2
通过3次及以下次数的谐波
另外记录图形，并标图号
滤除5次及以下次数的谐波
通过3次到5次的谐波
37 思考题
1．在实验箱采样频率fs固定为10KHz的条件下，要观察方波信号频带宽度内的各个谐波分量，方波信号的频率最高不能超过多少，为什么？
2．硬件实验内容7中输出信号各个谐波分量，与原来方波信号同样谐波分量相比，有没有发生失真？主要发生了什么类型的失真？为什么？
4 窗函数法FIR滤波器设计实验
41 实验目的
1．通过实验加深对FIR滤波器基本原理的理解。
2．学习使用窗函数法设计FIR滤波器，了解窗函数的形式和长度对滤波器性能的影响。
42 实验仪器
1．YBLD智能综合信号源测试仪 1台
2．双踪示波器 1台
3．MCOM－TG305数字信号处理与现代通信技术实验箱 1台
4．PC机（装有MATLAB、MCOM－TG305配套实验软件） 1台
43 实验原理
数字滤波器的设计是数字信号处理中的一个重要内容。数字滤波器设计包括FIR（有限单位脉冲响应）滤波器与IIR（无限单位脉冲响应）滤波器两种。
与IIR滤波器相比，FIR滤波器在保证幅度特性满足技术要求的同时，很容易做到严格的线性相位特性。设FIR滤波器单位脉冲响应h(n)长度为N，其系统函数H（z）为：
H（z）是z－1的N－1次多项式，它在z平面上有N－1个零点，原点z=0是N－1阶重极点，因此H（z）是永远稳定的。稳定和线性相位特性是FIR滤波器突出的优点。
FIR滤波器的设计任务是选择有限长度的h(n)。使传输函数H( )满足技术要求。FIR滤波器的设计方法有多种，如窗函数法、频率采样法及其它各种优化设计方法，本实验介绍窗函数法的FIR滤波器设计。
窗函数法是使用矩形窗、三角窗、巴特利特窗、汉明窗、汉宁窗和布莱克曼窗等设计出标准响应的高通、低通、带通和带阻FIR滤波器。
一、firl函数的使用
在MATLAB下设计标准响应FIR滤波器可使用firl函数。firl函数以经典方法实现加窗线性相位FIR滤波器设计，它可以设计出标准的低通、带通、高通和带阻滤波器。firl函数的用法为：
b=firl(n,Wn,/ftype/,Window)
各个参数的含义如下：
b—滤波器系数。对于一个n阶的FIR滤波器，其n+1个滤波器系数可表示为：b(z)=b(1)+b(2)z－1+…+b(n+1)z－n。
n—滤波器阶数。
Wn—截止频率，0≤Wn≤1，Wn=1对应于采样频率的一半。当设计带通和带阻滤波器时，Wn=[W1 W2],W1≤ω≤W2。
ftype—当指定ftype时，可设计高通和带阻滤波器。Ftype=high时，设计高通FIR滤波器；ftype=stop时设计带阻FIR滤波器。低通和带通FIR滤波器无需输入ftype参数。
Window—窗函数。窗函数的长度应等于FIR滤波器系数个数，即阶数n+1。
二、窗函数的使用
在MATLAB下，这些窗函数分别为：
1．矩形窗：w=boxcar(n)，产生一个n点的矩形窗函数。
2．三角窗：w=triang(n)，产生一个n点的三角窗函数。
当n为奇数时，三角窗系数为w(k)=
当n为偶数时，三角窗系数为w(k)=
3．巴特利特窗：w=Bartlett(n)，产生一个n点的巴特利特窗函数。
巴特利特窗系数为w(k)=
巴特利特窗与三角窗非常相似。巴特利特窗在取样点1和n上总以零结束，而三角窗在这些点上并不为零。实际上，当n为奇数时bartlett(n)的中心n－2个点等效于triang(n－2)。
4．汉明窗：w=hamming(n)，产生一个n点的汉明窗函数。
汉明窗系数为w(k+1)=054－046cos（ ） k=0,…,n－1
5．汉宁窗：w=hanning(n)，产生一个n点的汉宁窗函数。
汉宁窗系数为w(k)=05[1－cos( )] k=1,…,n
6．布莱克曼窗：w=Blackman(n),产生一个n点的布莱克曼窗函数。
布莱克曼窗系数为w(k)=042－05cos(2π )+08cos（4π ）] k=1,…,n
与等长度的汉明窗和汉宁窗相比，布莱克曼窗的主瓣稍宽，旁瓣稍低。
7．凯泽窗：w=Kaiser(n,beta)，产生一个n点的凯泽窗数，其中beta为影响窗函数旁瓣的β参数，其最小的旁瓣抑制α与β的关系为:
01102(α－087) α>50
β= 05842(α－21)04+007886(α－21) 21≤α≤50
0 α<21
增加β可使主瓣变宽，旁瓣的幅度降低。
8．契比雪夫窗：w=chebwin(n,r)产生一个n点的契比雪夫窗函数。其傅里叶变换后的旁瓣波纹低于主瓣r个db数。
44 实验内容
1．软件仿真实验：编写并调试MATLAB程序，观察不同窗，不同类型滤波器不同点数等共4种FIR滤波器的h(n)，并记录幅频特性和相频特性。
2．硬件实验：用窗函数法设计标准响应的FIR滤波器，在计算机上观察窗函数幅频特性、幅频特性和相频特性，然后下载到实验箱。用示波器观察输入输出波形，测试滤波器的幅频响应特性。
45 MATLAB参考程序和仿真内容
%%
%mode: 模式（1--高通；2--低通；3--带通；4--带阻）
%n: 阶数，加窗的点数为阶数加1
%fp: 高通和低通时指示截止频率，带通和带阻时指示下限频率
%fs: 带通和带阻时指示上限频率
%window:加窗（1--矩形窗；2--三角窗；3--巴特利特窗；4--汉明窗；
% 5--汉宁窗；6--布莱克曼窗；7--凯泽窗；8--契比雪夫窗）
%r: 代表加chebyshev窗的r值和加kaiser窗时的beta值
%sample: 采样率
%h: 返回设计好的FIR滤波器系数
%%
%mode: 模式（1--高通；2--低通；3--带通；4--带阻）
%n: 阶数，加窗的点数为阶数加1
%fp: 高通和低通时指示截止频率，带通和带阻时指示下限频率
%fs:
在公式（18）中，和分别表示状态变量与其修正量组成的列向量；为方阵，一般叫作雅可比矩阵，第i行j列元素为 ，它的大小为第i个函数对第j个变量求偏导；k则表示阵元素都在处取；同时，F(X)是由n个函数组成的n维列向量；在极坐标下，节点电压可如下表示：
在这里插入描述
（19）
若和为已知大小的功率，与从节点电压求得的有功和无功功率之差，为功率的不平衡量，则节点功率不平衡量可用如下公式计算：
在这里插入描述
（20）
节点功率可用各节点电压模值与相位表示，如下公式所示：
在这里插入描述
（21）
式（21）中，为节点i和j的相位差。
由以公式（18）-（21）推得牛顿法下，其潮流计算方程可写为：
在这里插入描述
（22）
公式（22）中，雅可比矩阵的各元素为
在这里插入描述
（23）
（24）
（25）
（26）
在这里插入描述
（27）
（28）
（29）
（30）
其中，节点导纳矩阵的元素由Gij 、Bij表示。
随着国内外配电系统自动化水平不断提高，电力行业人员也开始更加深入地研究配电网系统。配电网潮流计算作为DMS（配电管理系统）的重要基础，受到广大行业界人士的关注。因此，配电网潮流计算，已然成为配电网分析的重要内容。配电网与输电网相比，两者有明显不同，前者一般采用网格结构，线路参数R/X的值较大，三相负荷不对称程度明显。这些特点使得在输电网中计算有效，如牛顿法，不再适用于配电网。为此，有学者提出了适用于配电网的潮流算法，主要包括基于回路方程的潮流算法、前推回推法和改进的牛顿-拉夫逊法[17]（简称改进的牛拉法）。其中，基于回路方程的方法具有较强的网格处理能力和良好的收敛性，但该方法的节点数和分支数处理非常复杂。前推回推法是针对配电网的树状特性，可以避免潮流计算中的病态条件，同时速度更快。然而，由于其公式和算法与牛顿潮流算法不同，其在其它方面（如潮流优化）的应用将受到限制。
改进牛顿法通过对传统法进行一定的近似，将J阵写成UDUT 的形式。U仅由网络拓扑决定，是一个上三角矩阵；D是一个对角矩阵。在牛拉法中，需要对J阵因子分解与前代回代，改进法则只有前推回代的计算过程。它很好地改善了传统法以及前推回推法。经过算例计算结果证明，改进法可以避免J阵病态，且拥有前推回代法的收敛速度、精度，又由于它属于牛顿型算法，所以该算法已经得到了广泛的运用[18]。
下面附带电力系统分析牛顿法算例及matlab程序：
网络结构如下：系统结构图
系统参数如下：
在上图所示的简单电力系统中，系统中节点1、2为PQ节点，节点3为PV节点，节点4为平衡节点，已给定P1s+jQ1s=-030-j018 P2s+jQ2s=-055-j013 P3s=05 V3s=110 V4s=105∠0°
容许误差ε=10-5
节点导纳矩阵：
导纳矩阵
各节点电压：
节点 e f v ζ
10984637 -0008596 0984675 -0500172
20958690 -0108387 0964798 -6450306
31092415 0128955 1100000 6732347
41050000 0000000 1050000 0000000
各节点功率：
节点 P Q
1-0300000 -0180000
2–0550000 -0130000
3 0500000 -0551305
4 0367883 0264698
matlab程序如下：
// 牛顿法潮流计算matlab程序
clc;
Y=[1042093-8242876i -0588235+2352941i 3666667i -0453858+1891074i;
-0588235+2352941i 1069005-4727377i 0 -0480769+2403846i;
3666667i 0 -3333333i 0;
-0453858+1891074i -0480769+2403846i 0 0934627-4261590i];
%导纳矩阵
e=[1 1 11 105];%初始电压
f=zeros(4,1);
V=zeros(4,1);%节点电压
Ws=[-03 ; -018 ; -055 ; -013 ; 05 ; 11];%初始功率
W=zeros(6,1);
n=length(Y);%节点数
J=zeros(2(n-1));%雅可比矩阵
delta_v=zeros(1,6);
delta_w=Ws;
G=real(Y);
B=imag(Y);
S=zeros(4,2);
c=0;%循环次数
m=input('请输入PQ节点数：');
while max(abs(delta_w))>10^-5
for i=1:(n-1)%以下为求取雅可比矩阵
for j=1:(n-1)
if (i~=j)
J(2i-1,2j-1)=-(G(i,j)e(i)+B(i,j)f(i));
J(2i,2j)=-J(2i-1,2j-1);
J(2i-1,2j)=B(i,j)e(i)-G(i,j)f(i);
J(2i,2j-1)=J(2i-1,2j);
end
end
end
for j=1:(n-2)
J(6,2j-1)=0;
J(6,2j)=0;
end%以上为非对角线元素
s1=0;
s2=0;
for i=1:(n-1)
for j=1:n
s1=s1+(G(i,j)e(j)-B(i,j)f(j));
s2=s2+(G(i,j)f(j)+B(i,j)e(j));
end
J(2i-1,2i-1)=-s1-G(i,i) e(i)-B(i,i)f(i);
J(2i-1,2i)=-s2+B(i,i) e(i)-G(i,i)f(i);
s1=0;
s2=0;
end
for i=1:m
for j=1:n
s1=s1+G(i,j)f(j)+B(i,j)e(j);
s2=s2+(G(i,j)e(j)-B(i,j)f(j));
end
J(2i,2i-1)=s1+B(i,i) e(i)-G(i,i)f(i);
J(2i,2i)=-s2+G(i,i) e(i)+B(i,i)f(i);
s1=0;
s2=0;
end
J(6,5)=-2e(3);
J(6,6)=-2f(3);%对角线元素求解
for i=1:m
for j=1:n
s1=s1+e(i)(G(i,j)e(j)-B(i,j)f(j))+f(i)(G(i,j)f(j)+B(i,j)e(j));
s2=s2+f(i)(G(i,j)e(j)-B(i,j)f(j))-e(i)(G(i,j)f(j)+B(i,j)e(j));
end
delta_w(2i-1)=Ws(2i-1)-s1;
delta_w(2i)=Ws(2i)-s2;
W(2i-1)=s1;
W(2i)=s2;
s1=0;
s2=0;
end
for j=1:n
s1=s1+e(3)(G(3,j)e(j)-B(3,j)f(j))+f(3)(G(3,j)f(j)+B(3,j)e(j));
end
delta_w(5)=Ws(5)-s1;
delta_w(6)=(Ws(6)^2-(e(3)^2+f(3)^2));
W(5)=s1;
W(6)=sqrt(e(3)^2+f(3)^2);%以上求功率差值
delta_v=-inv(J)delta_w;
for i=1:(n-1)
e(i)=e(i)+delta_v(2i-1);
f(i)=f(i)+delta_v(2i);
end%求电压差值
c=c+1;
end
for x=1:4
V(x)=e(x)+f(x)1i;
end%节点电压
s1=0;
for x=3:4
for j=1:4
s1=s1+conj(Y(x,j))conj(V(j));
end
S(x,1)=real(V(x)s1);
S(x,2)=imag(V(x)s1);
s1=0;
end%PV与平衡节点功率
for x=1:2
S(x,1)=W(2x-1);
S(x,2)=W(2x);
end%节点功率
c
J
V
S
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
运行结果如下：
潮流计算结果
文章知识点与官方知识档案匹配
算法技能树首页概览
34356 人正在系统学习中
电磁MATLAB
微信公众号
程序从vx公众号获取，助力科研！
打开CSDN APP，看更多技术内容
matlab潮流程序,潮流计算的MATLAB源程序_weixin_39760433的博客-CSDN
1、简单潮流计算的牛顿拉夫逊程序,相关的原始数据数据数据输入格式如下:%B1是支路参数矩阵,第一列和第二列是节点编号。节点编号由小到大编写%对于含有变压器的支路,第一列为低压侧节点编号,第二列为高压侧节点%编号,将变压器的串联阻抗
继续访问
编程实现_某只旺仔的博客_电力系统潮流计算程序matlab
潮流计算是电力系统中最基本,应用最广泛的一种计算,是电力系统稳定计算和故障分析的基础。这个代码通过matlab实现了用极坐标表示的牛顿法和P-Q分解法进行潮流计算,使用IEEE14节点系统进行测试,计算结果和应用matpower的潮流计算完全一致。代码
继续访问
<em>matlab</em>导数计算
<p>matlab导数计算，通过编程实战掌握具体应用。包括matlab导数计算前、matlab导数计算中、matlab导数计算后。</p>
继续访问
(完整版)基于MATLAB牛顿拉夫逊法进行潮流计算doc
> %本程序的功能是用牛顿拉夫逊法进行潮流计算 n=input' 请输入节点数 :n=; nl=input' 请输入支路数 :nl=; isb=input' 请输入平衡母线节点号 :isb=; pr=input' 请输入误差精度 :pr=; B1=input' 请输入由各支路参数形成的矩阵 :B1=; B2=input' 请输入各节点参数形成的矩阵 :B2=; Y=zeros(n; e=zeros
牛顿拉夫逊潮流计算matlab程序
基于牛顿拉夫逊进行潮流计算，求得各节点电压，各支路功率流动，内附程序输入说明以及案例。
珍藏多年的matlab潮流计算程序源代码集合，包含多个潮流计算程序
达摩老生出品，必属精品，亲测校正，质量保证 资源名：珍藏多年的matlab潮流计算程序源代码集合，包含多个潮流计算程序 资源类型：matlab项目全套源码 源码说明： 全部项目源码都是经过测试校正后百分百成功运行的，如果您下载后不能运行可联系我进行指导或者更换。 适合人群：新手及有一定经验的开发人员
牛顿-拉夫逊法潮流计算——采用极坐标形式
采用极坐标形式的Newton-Raphson法进行潮流计算，并提供IEEE 14、57节点及New England 39节点的计算数据。使用方法：运行/修改PowerFlow_Newton_Polarm文件。
电力系统牛拉法潮流计算MATLAB程序
实现潮流计算的MATLAB代码，使用牛顿拉夫逊法，可更改结点数，支路数和导纳阵，通用性强
课设/毕业设计电力系统潮流计算（Matlab代码实现）
潮流计算是电力系统分析中的一种最基本的计算，对给定系统进行潮流计算可以得到各母线上的电压、网络中的功率分布及功率损耗等。本文介绍了潮流计算在电力系统分析中的作用及潮流计算的发展状况。通过对潮流计算所用的数学模型进行分析，建立潮流计算的基本方程。牛顿-拉夫逊法是目前广泛应用的一种潮流计算方法，本文阐述了牛顿-拉夫逊潮流计算的基本原理。详细展示了利用matlab平台编写潮流计算程序的具体过程，通过实例证明基于牛顿-拉夫逊法的潮流计算程序具有收敛速度快、占用内存小的优点。潮流计算的程序实现手段。
继续访问
传统牛顿拉夫逊法潮流计算matlab程序zip
采用何仰赞《电力系统分析》中介绍的牛顿拉夫逊法matlab编程，程序还能实现N-1校核和线路网损分析，适合学习电气工程领域的学子。
基于matpower的电力系统潮流计算matlab程序
它是由美国康奈尔大学电力系统工程研究中心（PSERC of Cornell University）的RAY D Zimmerman、Carlos E Murillo-Sánchez和甘德强在Robert J Thomas的指导下开发出来的，本章介绍的是MATPOWER40。每一个电网用变量名为“mpc”的结构体（structures）来定义，结构体mpc的不同字段用baseMVA、bus、branch、gen等来定义和返回电网的具体参数。列的数据类似于标准的IEEE 和PTI 列的数据格式。
继续访问
直流潮流计算matlab程序
直流潮流发的特点是用电力系统的交流潮流（有功功率和无功功率）等值的直流电流来代替。甚至只用直流电路的解析法来分析电力系统的有功潮流，而不考虑无功分布对有功的影响。这样一来计算速度加快，但计算的准确度有所降低，本方法适用于对潮流计算准确度要求不高的计算场景。θ为网络中各节点的电压相位角的向量；P为节点注入的有功功率向量。这就相当于线路两端的直流电位分别为θi和θj。（2）按照标幺值计算时，节点电压与其额定电压相差不大，故有：Ui≈Uj≈10；以IEEE9节点系统为算例，系统参数如下。
继续访问
牛顿-拉夫逊法潮流计算matlab程序,牛顿—拉夫逊法潮流计算MATLAB程序doc
！！！！！！！！！！！！########################牛顿—拉夫逊法潮流计算程序By Yuluo%牛顿--拉夫逊法进行潮流计算n=input('请输入节点数：n=');n1=input('请输入支路数：n1=');isb=input('请输入平衡母线节点号：isb=');pr=input('请输入误差精度：pr=');B1=input('请输入由支路参数形成的矩阵：B1=');B
继续访问
matlab程序 潮流计算,潮流计算matlab程序
《潮流计算matlab程序》由会员分享，可在线阅读，更多相关《潮流计算matlab程序(3页珍藏版)》请在人人文库网上搜索。1、clear;%各节点参数：节点编号，类型，电压幅值，电压相位，注入有功，注入无功 %类型：1PQ节点，2PV节点，3平衡节点%本程序中将最后一个节点设为平衡节点R_1=1 1 10 0 02 02j;2 1 10 0 -045 -015j;3 1 10 0
继续访问
matlab潮流计算函数,基于牛顿拉夫逊法潮流计算的matlab实验报告(含源程序和结果)
《基于牛顿拉夫逊法潮流计算的matlab实验报告(含源程序和结果)》由会员分享，可在线阅读，更多相关《基于牛顿拉夫逊法潮流计算的matlab实验报告(含源程序和结果)(24页珍藏版)》请在人人文库网上搜索。1、基于牛顿拉夫逊法潮流计算的matlab实验报告一、 实验目的和要求1学习掌握matlab的基本用法2应用MATLAB语言编写具有一定通用性的牛顿-拉夫逊法潮流计算程序。要求：(1)潮流计
继续访问
matlab的做潮流计算,Matlab实现潮流计算程序
关于Matlab趋势计算programdoc文件的实现，爱文共享的信息具有丰富的相关文档，每天都有成千上万的行业名人在该站点上共享最新信息程序代码如下: 读取数据clcclearfilename ='txt'a = textread(filename)n = a(，)pinghengjd = a(，)phjddianya = a(，)jingdu = a(，)b = zeros (，)j =
继续访问
天然气潮流计算matlab程序
此外，由于天然气在管道内传输时会产生压力损耗，所以每隔一段都要设置压缩机来提升节点压力，保证天然气系统的正常运行。（3）如果节点中存在压缩机，则可以求得压缩机消耗的等效天然气流量，将其作为一个负荷加入到潮流计算中；管道或者节点的约束，例如流量、节点压力、压缩机消耗电能和压缩比等，如下面四式所示。（2）根据天然气节点压力与节点间流量的关系式求出前一个节点的节点压力；（1）由用户负荷求得微型燃气轮机的功率
欢迎分享，转载请注明来源：内存溢出
微信扫一扫
支付宝扫一扫
评论列表（0条）