福建工程学院电网电流采集与处理.docx
《福建工程学院电网电流采集与处理.docx》由会员分享，可在线阅读，更多相关《福建工程学院电网电流采集与处理.docx（35页珍藏版）》请在冰豆网上搜索。
福建工程学院电网电流采集与处理
第1章课程设计任务书
一、课程设计的目的
（1）进一步理解数字信号处理的基本概念、基本理论和基本方法。
（2）熟悉在MATLAB环境下电网负载电流采集的方法。
（3）学会用MATLAB软件对信号进行分析和处理。
（4）综合运用数字信号处理理论知识，掌握用MATLAB软件设计FIR和IIR数字滤波器的方法。
（5）提高依据所学知识及查阅的课外资料来分析问题解决问题的能力。
二、课程设计的要求
（1）认真阅读课程设计任务书，熟悉有关设计资料及参考资料，熟悉各种设计规范的有关内容，认真完成任务书规定的设计内容。
（2）学生在规定的时间内独立完成规定的内容和工作量。
（3）课程设计的成果为“课程设计任务书”。
要求计算准确、文字通顺、图面布置合理、正确清晰、符合制图标准及有关规定。
（4）课程设计任务书是一项工程或事件的真实情况的反映，要用事实、数据说话。
任务书不需要目录、摘要、关键词、脚注、参考文献等。
课程设计任务书应在2000字以上（含相关图纸和计划书等），打印成稿。
三、设计的内容
1、电网负载电流采集
三相可控整流器带纯电阻负载（负载可变），相控角可设置，观察波形，并记录，仿真时间0.2秒，记录的数组为iLoad。
2、负载电流频谱分析
在Matlab中，利用fft函数对信号进行频谱分析，得到信号的频谱特性。
3、设计数字滤波器
给定滤波器的性能指标如下：
低通：
fp=60Hz，fs=80Hz，Ap=1dB，As=45dB。
高通：
fs=370Hz，fp=390Hz，Ap=1dB，As=45dB。
带通：
fs1=210Hz，fp1=230Hz，fp2=270Hz，fs2=290Hz，Ap=1dB，As=45dB。
采用窗函数法和双线性变换法设计上面要求的3种滤波器，并画出滤波器的频率响应，并根据幅频特性曲线对所设计的数字滤波器进行可行性分析。
4、用滤波器对信号进行滤波
用自己设计的滤波器对采集到的信号进行滤波处理，画出滤波前后信号的时域波形及频谱，并对滤波前后的信号进行对比，分析信号的变化。
5、在CCS构造3个频率信号叠加，f1=50Hz，f2=50*5Hz，f3=50*7Hz，编制软件，实现FIR低通和带通滤波，低通滤波器要滤除两个高频信号，带通滤波器要滤除低频信号和高频信号，Ap=1dB，As=50dB。
用CCS画图工具画出滤波前及滤波后信号波形。
（先设计FIR滤波器，然后对滤波器系数定标并移植）
第2章设计原理
一、窗函数法设计原理
希望逼近的滤波器频率响应函数为：
单位脉冲响应为：
hd（n）是中心点在α的偶对称无限长非因果序列，如下图（a）。
要得到有限长的h（n），必须加窗截断，矩形窗RN（n）如下图（b）。
按照第一类线性相位FIR滤波器的约束条件，h（n）必须偶对称，对称中心应为（N-1）/2，因此，必须取α=（N-1）/2，所以有
窗函数设计法的时域波形（矩形窗，N=30）
加窗截断产生的影响
实际滤波器的幅度函数为
对实际FIR滤波器频率响应的幅度函数起影响的是窗函数频率响应的幅度函数
用一个有限长序列h（n）代替hd（n），会引起截断效应，体现在过渡带加宽以及通带和阻带内的波动，阻带衰减减小。
二、用窗函数法设计FIR数字滤波器的步骤
（1）根据阻带最小衰减选择窗函数类型，根据过渡带的指标要求估计窗口长度Ｎ
（2）构造希望逼近的频率响应函数
（3）计算hd（n）
（4）加窗
，得到计算结果
（5）求
，检验是否满足设计要求，如果不满足，重新设计
三、窗函数法的MATLAB设计函数
b=fir1（M,wc,’ftype’,window）
M：
FIR滤波器阶数，M=N-1
window=kaiser（N,alph），凯塞窗
wc:
对pi归一化的数字频率
wc=（wp+ws）/2，ftype不填，设计FIR低通滤波器
wc=（wp+ws）/2，ftype=high，设计FIR高通滤波器
wc=[wcl,wcu]，ftype不填，设计FIR带通滤波器
四、用双线性变换法设计IIR数字滤波器
无限脉冲响应数字滤波器
[N,wc]=buttord（wp,ws,Rp,As,’s’）
计算Butterworth模拟滤波器的阶数N和3dB截止频率wc。
wp，ws和wc是实际的模拟角频率（单位：
rad/s）。
[B,A]=butter（N,wc,'ftype','s'）
计算Butterworth模拟滤波器系统的分子和分母多项式的系数向量B和A。
由系数向量B和A写出模拟滤波器的系统函数H（s）为：
[N,wpo]=ellipord（wp,ws,Rp,As,‘s’）%wpo为通带边界频率
[B,A]=ellip（N,Rp,As,wpo,'ftype','s'）
模拟高通、带通、带阻滤波器的设计波形图
第3章设计内容
（1）电网负载电流采集
1.总的连接电路图如下：
2.三相电网电路图
3.三相可控整流器带纯电阻负载电路图如下：
（2）观察波形
1.电流采集结果iLoad波形
2.iL波形
3.IL1波形
4.设置全局变量负载电流一个电网周期采样点数N和采样周期Ts
（3）负载电流频谱分析
1.程序如下：
%课程设计---获取离散信号并做频谱分析，了解对象的频率分布情况
closeall;
clear;
clc;
loadiLoad;%读取采样结果
Fs=6.4*1000;%采样频率6.4kHz
Ts=1/Fs;%采样周期
N=0.02/Ts;%负载电流一个电网周期采样点数
OrdMax=22;%阶数最大值
t=iLoad.time;%负载电流时间
L=length（t）%负载电流总采样点数
current=iLoad.signals.values;%负载电流等间隔采样，得到离散信号:
current
（1）,current
（2）,current（3）,...,current（L）
xn=current（N:
-1:
1）;%一个电网周期的样本
Xk=fft（xn）;
%以下是绘图部分
figure
（1）
subplot（3,1,1）,plot（current）;title（'负载电流'）;
xlabel（'t/s'）;ylabel（'i/A'）;%物理量及单位
gridon;%加网格
xlim（[0,L]）;ylim（[-60,60]）;%坐标限定
subplot（3,1,2）,stem（abs（Xk）/N*2,'.'）;title（'负载电流频谱'）;%plot（abs（Xk）/N*2）;
gridon;
xlim（[0,N]）;ylim（[0,60]）;
k=0:
22;
subplot（3,1,3）,stem（k,abs（Xk（k+1））/N*2,'.'）;
title（'负载电流频谱'）;%plot（k,abs（Xk（k+1））/N*2）;
gridon;
xlim（[0,22]）;ylim（[0,60]）;
2.利用fft函数对信号进行频谱分析，得到信号的频谱特性如下
（4）根据给定滤波器的性能指标设计FIR低通、高通、带通滤波器
并用滤波器对信号进行滤波
1.程序如下：
%**************************加载FIR滤波器*************************%
%*****************************低通滤波器**************************%
fp=60;fs=80;rs=21;%设计指标
wp=（2*pi*fp）/Fs;
ws=（2*pi*fs）/Fs;
Bt=ws-wp;
alph=0.5842*（rs-21）^0.4+0.07886*（rs-21）;%计算kaiser窗的形状参数
M=ceil（（rs-8）/2.285/Bt）;%根据（7.2.17）式计算kaiser窗所需阶数M
wc=（wp+ws）/2/pi;%计算理想低通滤波器通带截止频率（关于π归一化），即wc:
对pi归一化的数字频率
h1=fir1（M,wc,kaiser（M+1,alph））;%调用kaiser计算低通FIRDF的h（n）
hh1=round（h1*32768）;
fld=fopen（'hh1.txt','wt'）;
fprintf（fld,'%d,',hh1）;
fclose（fld）;
%以下是绘图部分
figure
（2）
n=1:
M+1;
subplot（3,2,1）;stem（n,h1,'.'）;gridon;
xlim（[0,M+1]）;%ylim（[-0.05,0.05]）;
title（'FIR低通滤波器单位脉冲响应h（n）'）;xlabel（'n'）;ylabel（'h（n）'）;
N_DFT=512;
[H,w]=freqz（h1,1,N_DFT）;%求幅频响应
k=1:
N_DFT;
subplot（3,2,2）;plot（w/pi,20*log10（abs（H（k））））;title（'幅频响应'）;
gridon;xlabel（'\omega/\pi'）;ylabel（'Magnitude（dB）'）;
%3.2.1低通滤波器应用（对负载电流进行滤波及滤波后的频谱分析）
current1=fftfilt（h1,current）;%对负载电流进行滤波
fft1=fft（current1（2*N+1:
2*N+1+N））;%从第二个周期开始取一个周期滤波后的电流值
n=1:
L;k=0:
22;
figure
（2）
subplot（3,2,3）;plot（t,current（n））;gridon;
ylim（[-70,70]）;
title（'原始电流信号'）;xlabel（'t/s'）;ylabel（'i/A'）;
subplot（3,2,5）;plot（t,current1（n））;gridon;
xlim（[0,0.2]）;ylim（[-70,70]）;
title（'FIR低通滤波后电流信号'）;xlabel（'t/s'）;ylabel（'i/A'）;
subplot（3,2,4）;stem（k,abs（Xk（k+1））/N*2,'.'）;gridon;%plot（k,abs（Xk（k+1））/N*2）;
xlim（[0,22]）;ylim（[0,60]）;
title（'原始电流信号频谱'）;xlabel（'k'）;ylabel（'谐波幅值'）;
subplot（3,2,6）;stem（k,abs（fft1（k+1））/N*2,'.'）;gridon;%plot（k,abs（fft1（k+1））/N*2）;
xlim（[0,22]）;ylim（[0,70]）;
title（'FIR低通滤波后电流信号频谱'）;xlabel（'k'）;ylabel（'谐波幅值'）;
%*****************************高通滤波器**************************%
%3.2高通滤波器设计（Kaiser窗）
fs=370;fp=390;rp=1;rs=21;%设计指标----8次以下谐波滤除
ws=2*pi*fs/Fs;wp=2*pi*fp/Fs;
Bt=wp-ws;
%alph=0.112*（rs-8.7）;%根据（7.2.16）式计算kaiser窗的参数α
alph=0.5842*（rs-21）^0.4+0.07866*（rs-21）;
M=ceil（（rs-7.95）/2.286/Bt）+1;%根据（7.2.17）式计算kaiser窗所需阶数M
wc=（wp+ws）/2/pi;%计算理想高通滤波器通带截止频率（对π归一化）
h2=fir1（M,wc,'high',kaiser（M+1,alph））;%所设计的FIR滤波器---单位脉冲响应h（n）
hh2=round（h2*32768）;
fld=fopen（'hh2.txt','wt'）;
fprintf（fld,'%d,',hh2）;
fclose（fld）;
%以下是绘图部分
figure（3）
n=1:
M+1;
subplot（3,2,1）;stem（n,h2,'.'）;gridon;
xlim（[0,M+1]）;%ylim（[-0.05,0.05]）;
title（'FIR高通滤波器单位脉冲响应h（n）'）;xlabel（'n'）;ylabel（'h（n）'）;
N_DFT=512;
[H,w]=freqz（h2,1,N_DFT）;
k=1:
N_DFT;
subplot（3,2,2）;plot（w/pi,20*log10（abs（H（k））））;title（'幅频响应'）;
gridon;xlabel（'\omega/\pi'）;ylabel（'Magnitude（dB）'）;
%3.2.1高通滤波器应用（对负载电流进行滤波及滤波后的频谱分析）
current2=fftfilt（h2,current）;
fft2=fft（current2（2*N+1:
2*N+1+N））;
n=1:
L;k=0:
22;
figure（3）
subplot（3,2,3）;plot（t,current（n））;gridon;
ylim（[-70,70]）;
title（'原始电流信号'）;xlabel（'t/s'）;ylabel（'i/A'）;
subplot（3,2,5）;plot（t,current2（n））;gridon;
xlim（[0,0.2]）;ylim（[-70,70]）;
title（'FIR高通滤波后电流信号'）;xlabel（'t/s'）;ylabel（'i/A'）;
subplot（3,2,4）;stem（k,abs（Xk（k+1））/N*2,'.'）;gridon;%plot（k,abs（Xk（k+1））/N*2）;
xlim（[0,22]）;ylim（[0,60]）;
title（'原始电流信号频谱'）;xlabel（'k'）;ylabel（'谐波幅值'）;
subplot（3,2,6）;stem（k,abs（fft2（k+1））/N*2,'.'）;gridon;%plot（k,abs（fft2（k+1））/N*2）;
xlim（[0,22]）;ylim（[0,70]）;
title（'FIR高通滤波后电流信号频谱'）;xlabel（'k'）;ylabel（'谐波幅值'）;
%*****************************带通滤波器**************************%
%3.2带通滤波器设计（Kaiser窗）
fls=210;flp=230;%设计指标
fup=270;fus=290;
rs=21;
wlp=2*pi*flp/Fs;
wls=2*pi*fls/Fs;
wup=2*pi*fup/Fs;
wus=2*pi*fus/Fs;
Bt=wlp-wls;
alph=0.5842*（rs-21）^0.4+0.07886*（rs-21）;%计算kaiser窗的形状参数P210
M=ceil（（rs-8）/2.285/Bt）;
wc=[（wls+wlp）/2/pi,（wus+wup）/2/pi];
h3=fir1（M,wc,kaiser（M+1,alph））;
hh3=round（h3*32768）;
fld=fopen（'hh3.txt','wt'）;
fprintf（fld,'%d,',hh3）;
fclose（fld）;
%以下是绘图部分
figure（4）
n=1:
M+1;
subplot（3,2,1）;stem（n,h3,'.'）;gridon;
xlim（[0,M+1]）;%ylim（[-0.05,0.05]）;
title（'FIR带通滤波器单位脉冲响应h（n）'）;xlabel（'n'）;ylabel（'h（n）'）;
N_DFT=512;
[H,w]=freqz（h3,1,N_DFT）;
k=1:
N_DFT;
subplot（3,2,2）;plot（w/pi,20*log10（abs（H（k））））;title（'幅频响应'）;
gridon;xlabel（'\omega/\pi'）;ylabel（'Magnitude（dB）'）;
%3.2.1带通滤波器应用（对负载电流进行滤波及滤波后的频谱分析）
current3=fftfilt（h3,current）;
fft3=fft（current3（2*N+1:
2*N+1+N））;
n=1:
L;k=0:
22;
figure（4）
subplot（3,2,3）;plot（t,current（n））;gridon;
ylim（[-70,70]）;
title（'原始电流信号'）;xlabel（'t/s'）;ylabel（'i/A'）;
subplot（3,2,5）;plot（t,current3（n））;gridon;
xlim（[0,0.2]）;ylim（[-70,70]）;
title（'FIR带通滤波后电流信号'）;xlabel（'t/s'）;ylabel（'i/A'）;
subplot（3,2,4）;stem（k,abs（Xk（k+1））/N*2,'.'）;gridon;%plot（k,abs（Xk（k+1））/N*2）;
xlim（[0,22]）;ylim（[0,60]）;
title（'原始电流信号频谱'）;xlabel（'k'）;ylabel（'谐波幅值'）;
subplot（3,2,6）;stem（k,abs（fft3（k+1））/N*2,'.'）;gridon;%plot（k,abs（fft3（k+1））/N*2）;
xlim（[0,22]）;ylim（[0,70]）;
title（'FIR带通滤波后电流信号频谱'）;xlabel（'k'）;ylabel（'谐波幅值'）;
2.采用窗函数法设计上面要求的3种滤波器的频率响应
（a）低通滤波器
（b）高通滤波器
（c）带通滤波器
将幅频响应横坐标限定0.2，便于观察
（a）低通滤波器
（b）高通滤波器
（c）带通滤波器
3.根据幅频特性曲线对所设计的数字滤波器进行可行性分析
由幅频特性曲线可以验证各项指标满足要求，如不满足要求，则要改变N，或改变相关参数，然后重新计算。
4.滤波前后信号的时域波形及频谱如上图，并对滤波前后的信号进行对比，分析信号的变化。
我们知道图像中采集电流的谐波最高到50次，幅值差不多50几mA,主要能量集中在基波，谐波次数月光，幅值越小。
只研究有限项，把无限项过滤掉。
从三个滤波器的时域图像中只能看到单频信号，看不出那些频率分量。
从频域看低通保留1次谐波，高通滤除掉8次以下谐波，带通保留5次谐波。
（5）根据给定滤波器的性能指标设计低通、高通、带通滤波器并用滤波器对信号进行滤波
1.程序如下;
%**************************加载IIR滤波器*************************%
%***************************低通滤波器***************************%
%低通滤波器fp=60Hz，fs=80Hz，Ap=1dB，As=45dB。
fp=60;fs=80;rp=1;rs=45;%设计指标
wp=2*pi*fp;%数字角频率
ws=2*pi*fs;
[N1,wpo]=ellipord（wp,ws,rp,rs,'s'）;%计算低通模拟滤波器阶数N1和通带边界频率wpo
[B,A]=ellip（N1,rp,rs,wpo,'s'）;%计算低通模拟滤波器系统函数系数向量B，A
[Bz,Az]=bilinear（B,A,Fs）;%用双线性变换法转换成数字滤波器
%以下是绘图部分
figure（5）
subplot（3,2,1）;
fk=0:
（50*OrdMax）/512:
（50*OrdMax）;
wk=2*pi*fk;
Hk=freqs（B,A,wk）;
plot（fk,20*log10（abs（Hk）））;gridon;
title（'椭圆模拟低通滤波器幅频响应'）;xlabel（'频率（kHz）'）;ylabel（'幅度（db）'）;
xlim（[0,50*OrdMax]）;ylim（[-80,1]）;
subplot（3,2,2）;
wk=0:
pi/512:
pi;
Hz=freqz（Bz,Az,wk）;
plot（wk/pi,20*log10（abs（Hz）））;gridon;
title（'IIR椭圆数字低通滤波器幅频响应'）;xlabel（'\omega/\pi'）;ylabel（'幅度/db'）;
axis（[0,0.2,-80,1]）;
%低通滤波器应用（对负载电流进行滤波及滤波后的频谱分析）
current1=filter（Bz,Az,current）;
fft1=fft（current1（3*N+1:
3*N+1+N））;%从第三个周期开始取一个周期滤波后的电流值
n=1:
L;k=0:
22;
figure（5）
subplot（3,2,3）;plot（t,current（n））;gridon;
ylim（[-70,70]）;
title（'原始电流信号'）;xlabel（'t/s'）;ylabel（'i/A'）;
subplot（3,2,5）;plot（t,current1（n））;gridon;
xlim（[0,0.2]）;ylim（[-70,70]）;
title（'IIR低通滤波后电流信号'）;xlabel（'t/s'）;ylabel（'i/A'）;
subplot（3,2,4）;stem（k,abs（Xk（k+1））/N*2,'.'）;gridon;%plot（k,abs（Xk（k+1））/N*2）;
xlim（[0,22]）;ylim（[0,60]）;
title（'原始电流信号频谱'）;xlabel（'k'）;ylabel（'谐波幅值'）;
subplot（3,2,6）;stem（k,abs（fft1（k+1））/N*2,'.'）;gridon;%plot（k,abs（fft1（k+1））/N*2）;
xlim（[0,22]）;ylim（[0,60]）;
title（'IIR低通滤波后电流信号频谱'）;xlabel（'k'）;ylabel（'谐波幅值'）;
%****************************高通****************************%
fs=370;fp=390;rp=1;rs=45;
wp=2*pi*fp;%数字角频率
ws=2*pi*fs;
[N2,wpo]=ellipord（wp,ws,rp,rs,'s'）;
[B,A]=ellip（N2,rp,rs,wpo,'high','s'）;
[Bz,Az]=bilinear（B,A,Fs）;
%以
升级会员 