电网谐波测量

来源：爱够旅游网

1 绪论

随着国民经济的发展和人们生活水平的提高，电力电子产品广泛地应用于工业控制领域，用户对电能质量的要求也越来越高，其中最为突出的是电压质量和谐波的问题，因此，如何提高电压质量、治理谐波就成为输配电技术中最为迫切的问题之一。所以，面对我国目前电网结构薄弱和输配电技术普遍存在的技术手段的落后、自动化水平低的现状，针对电压质量和谐波问题，研究电网谐波治理问题和无功补偿新技术及新装备，具有十分重要的理论和现实意义[3]。

1．1 谐波的定义

“谐波”这一名词起源于声学，在声学中谐波表示一根弦或一个空气柱以基波频率的倍数频率振动。电气学中所谓电网谐波，就是电网正弦电压波形畸变后，其波形可以按傅立叶级数进行分解，除了基波（50HZ）之外，还有一系列频率为基波频率整数倍的正（余）弦波，这些正（余）弦波称之为谐波。正是由于这些谐波注入了电网，就使得电网电压波形畸变[14]。

1．2 谐波的危害

电网谐波的危害主要有以下几点：

1、相同频率的谐波电压余谐波电流要产生同此谐波的有功功率与无功功率，从而降低电网电压，浪费电网容量。

2、高次谐波能使电容器出现过电流与过负荷，温度增高，寿命减少，甚至出现发热、鼓肚、击穿或爆炸事故。同时在电压已经畸变的电网中，电容器的投入，还可能使电网的谐波加剧（谐波放大现象）。

3、谐波往往引起继电保护不工作或误动作，从而造成设备与系统的事故，尤其是半导体继电保护与整流型继电保护更为严重。

4、谐波能增大仪表的计量误差，干扰通讯网络的正常工作。

5、电机中有谐波电流，且频率接近某个零件的固有频率时，使电机产生机械振动并发出很大的噪声。

6、谐波对人体有影响。从人体生理学来看，人体细胞在受到刺激兴奋时，会在细胞膜静息电位基础上发生快速电波动或可逆翻转。其频率如果与谐波频率相接近，电网谐波的磁辐射就会直接影响人的脑磁场和心磁场。

1．3 谐波的产生

电网谐波来源于三个方面：其一是发电源质量不高产生谐波；其二是输电网产生谐波；其三是用电设备产生的谐波。其中以电气设备产生的谐波最多，具体情况如下：

1、整流设备。由于晶闸管整流的广泛应用（如电力机车的、路电解槽、电池充电器等），给电网造成大量的谐波。统计表明：由于整流装置产生的谐波占所有谐波的40%左右，这是最大的谐波源。

2、电弧炉、电石炉。由于加热原料时电炉的三项电极很难同时接触到高低不平的炉料，使得燃料不稳定，引起三项负荷不平衡，产生谐波电流，经变压器的△形连接线圈而注入电网。其中主要是2～7次的谐波，平均可达基波的8%～20%，最大可达45%。

3、电力变压器。由于变压器铁心的饱和，磁化曲线的非线性，加上设计变压器时考虑经济型，其工作磁密选择在磁化曲线的近饱和段上，这样就使得磁化电流呈尖顶波形，因而含有奇次谐波，其次谐波电流可达额定电流的0.5%。另外变压器空载合闸时出现的涵流中也含有大量的谐波量。

4、家用电器。如电视机、录像机、电子调光灯具、调温炊具等，因具有调压整流装置，会产生较深的奇次谐波；在洗衣机、电风扇、空调器等有绕组的设备中，因不平衡电流的变化也能是波形改变。

1．4 谐波抑制的三种方法

为了产少谐波的危害，根据谐波的来源，抑制谐波的方法大致从三个方面考虑： 1、确保电源的质量从发电机本身来说，产生谐波是很少的，这主要就是要从管理上保证发配电的质量，尤其是对于自发电设备的企业，要制定相应的措施予以保证。

2、减少输电网谐波的产生与磁辐射如尽可能采用高压直流输电；低压供电尽可能采用地下电缆；高层建筑内部的供电线应用铁管来防止磁辐射，还应做好大楼建筑的接地系统的设计与实施安装。

3、用电负荷产生的谐波前已述及，用电负荷是电网谐波的最主要的来源。首先要求用户对电网电压的畸变率与谐波电流的要做到负荷国家要求。根据国家水电部颁发的SD126-84规定,如果超过极限值，则应采取抑制措施。其次，就是要把维护电网不受谐波污染定成规章制度，用户与电气设备制造厂商都要配合发供电部门齐抓共管。具体来说应做到下面几点：

（1）对容量在100kVA以上的整流装置与非线性用电设备，必须就地增设滤波装置，吸收高次谐波。

（2）具有整流元件的设备，虽是小型用电器，也要增加整流相数与脉冲数，或把晶闸管移相触发改为过零触发，或者增设滤波装置，做到尽可能减少谐波。例如有些厂家生产的晶闸管调光台灯就增设了滤波原件。虽然成本有所增加，但能抑制谐波对电视机等设备的干扰，应受到推广。

（3）制造电磁感应式设备的厂家，要采用高质量的电磁材料，也可适当降低磁密。对于劣质元件的电磁设备要出售使用。

（4）对于电容器组等大容量的设备，应在回路内串电抗器或增设限流装置，抑制高次谐波。

（5）合理选择供电电压，尽量保持三相平衡。对能产生谐波的大容量设备还可采用较高电压供电，以减少系统的阻抗，从而减少谐波引起的网损。

（6）在产品的设计上还可采用干扰频率转嫁措施。如在变压器降压整流稳压的设备中，采用开关电源的逆变技术，就能大大降低低频段的污染，而转嫁到150kHz以上频率的污染，而这种较高频率的谐波容易被抑制。

总之，面对日趋严重的电网谐波污染，必须引起高度重视，这就需要发供电部门、用户及电气设备的制造厂商都从大局出发，共同努力，清扫和维护电网的“公共卫生”。同时只有这样，才能有利于各行各业的生产发展，才能有利于人民生活水平的不断提高。

2 课题分析

2．1 课题来源

随着电力电子技术的迅速发展，各种变频器、变流器、开关电源和电抗器的应用日益增多，电网随之产生的谐波污染也日益严重，影响了电力系统的安全运行及线路周围的电磁环境，因此电网谐波测量成为热门研究课题。

2．2 设计思路

将220V的交流电压通过电压采样环节，得到模拟量电压，将得到的模拟量电压接到USB2002采样卡模拟输入通道，接着由采样卡的采样环节变换成离散的数字量，然后运用傅立叶变换作谱分析，由于栅栏效应，谱线位置与分析点可能不重合，可采用具有实谱特性的窗函数对信号进行加窗处理，然后对信号及其时移连续作两次FFT分

析（实为一次复时域信号分析），则在窗谱的主瓣内将含有若干个由窗泄漏造成的相位相同的谱分析值，据此可以判定这些分析值由同一条真实谱线引起。继而可精确求出诸分立谱线的频率、幅度，然后将计算获得的数据显示在计算机屏幕上。

本课题是从电网中采集数据然后运用DFT来分析电网中存在的谐波。根据课题的需要可分成两大块，一块是数据采集，一块是数据处理。在数据采集模块里包括数据的采集和保存，本设计采用北京阿尔泰科贸有限公司生产的USB2002型号的采集卡进行采集。在数据处理模块中将采用DFT来进行数据处理从而得到谱线的幅度和频率，进而求出纹波系数。由于要用到编程来解决数据的采集、存盘和处理，我们将采用VC6.0做用户界面，通过调用MATLAB引擎的方法处理数据，以实现本课题的要求。

3 数据采集

3.1 数据采集卡的介绍

USB又称之为通用串行总线，它将计算机和外设连接在一起，能进行数字图象处理，同时为数字化设计提供了许多创造空间。

课题中我们所用的采集卡型号是USB2002。它是USB总线兼容的数据采集板，可经过USB电缆接入计算机，构成实验室、产品质量检验中心、野外测控、医疗设备等领域的数据采集、波形分析和处理系统，也可构成工业生产过程控制的系统。而且它具有体积小，即插即用等特点，因此是便捷式系统用户的最佳选择。现在我们使用的是USB2.0的版本，它的传输速度最高达到480MBit。同时本卡的使用是由软硬件共同完成的。

3.2 数据采集卡的技术指标及应用

1、模拟信号输入部分

模拟通道输入数：32路单端/16路双端模拟信号输入，支持通道扫描及伪同时采集（同步采集）扫描模式（注：伪同步模式：即模拟同步采集模式，采样被定时器或外部时钟启动后，USB2002以400KHZ的最大速度对用户设置的一组通道采样，结束后等待下一次启动，如此循环采样）。

（1）模拟电压输入范围：±5V，±10V（AD79-1）（2）模拟输入阻抗：100MΩ （3）模拟输入共模电压范围：>±2V

（4）放大器建立时间：2μS 2、A/D转换电路部分

（1）A/D分辨率：14Bit（16384）（2）非线性误差：±1LSB（最大）（3）转换时间：2.5μS

（4）系统测量精度（满量程）：0.1％（5）采样速度：400KHZ

3、应用

数据采集卡的主要应用在以下几个方面：野外采集、信号采集和医疗设备

3.3 程序调试

现在由于在实验室操作，为了安全问题，我们先将电网220V的电压通过变压器降压，然后再用电阻分压（当然使用了变压器会对我们要处理的电网谐波有一定的影响）。将电压采样电路和数据采集卡连好后，我们便将电源打开，接着打开我们要运行的程序。点击程序运行按钮，系统自动弹出一个框，在这个框中，我们进行一些参数的设置，参数设置完后我们便开始数据采集，采集的数据将在左半边出现，而图形以正弦波的形式在右半边出现。采集结束后点击存盘按钮，数据将存入相应的盘中。在操作中我们发现如果使用单通道输入会有一些干扰，而使用双通道输入则可以抑制一些干扰。

4 数据处理

4．1 综述

数据处理就是对采集到的信号进行分析，分析的步骤有以下几方面：首先取两组相差为一个采样周期的两组数据，然后对数据加布莱克曼－哈里斯窗，接着对数据作DFT变换，最后用查询方法寻找真实谱线，流程图如图4.1所示：

4．2 VC++调用MATLAB

4．2．1 VC调用MATLAB的三种方法

Visual C++是基于Windows平台下的一种应用程序开发环境,现在已升级到了6.0版本。应用VC6.0可以方便地开发各种用户需要的界面和工具,并且能够直接与系统的底层硬件交换数据。在VC6.0环境下开发的软件,具有容易维护、界面友好及运行速度

快等许多优点。但是, 在VC6.0下进行编程是比较烦琐的,不易掌握。因此, 将VC6.0

开始 VC调用MATLAB引擎取两组数据加窗作DFT变换查找真实谱线查找完？结束

图4.1 数据处理部分流程图

环境与MATLAB进行有效连接将会使得VC程序的开发大大简化[13]。

以下是三种VC调用MATLAB的方法，它们各具特色,在应用中要结合具体需求进行选择。第一种：MCC 编译器法可以直接对MATLAB 产生的.m文件进行编译,在MATLAB提供的C/C++函数库、图形库的支持下,能将绘制各种复杂图形的MATLAB程序高效的转化为C/C++代码,与单纯用VC编程实现相比大大减轻编程人员的工作量和编程难度；第二种：采用Matcom 编译法进行转换更加方便,且生成的C/C++代码的可读性较好,是解决图形、图像处理问题的一种理想方法；第三种：使用MATLAB计算引擎法时,由于其工作时MATLAB后台工作,只有小部分引擎通信函数库与程序相连,因此节省了资源,提

高了应用程序的整体性能和处理效率。

4．2．2 调用MATLAB引擎

MATLAB允许用户通过MATLAB引擎调用MATLAB的函数，即将MATLAB当作应用程序进行数据计算的引擎。MATLAB提供了一系列操作MATLAB引擎的API函数，通过这些API函数，用户可以充分发挥MATLAB进行矩阵计算的优势，将应用程序的计算功能交给MATLAB引擎完成，而界面部分则可以通过采用VC++来实现。这些为应用程序的实现提供了很大的灵活性。

1、API函数介绍

1）int engClose(Engine *ep); 退出MATLAB引擎。

2）int engEvalString(Engine *ep,const char *string); 使MATLAB引擎执行字符串string中的表达式。

3）mxArray *engGetVariable(Engine *ep, const char *name); 从MATLAB引擎工作空间中复制名字为name的变量。 4）int engGetVisible(Engine *ep,bool *value); 判断MATLAB引擎工作窗口是否可见。

5）Engine *engOpen(const char *startcmd);

启动一个MATLAB引擎，在Windows操作环境下startcmd参数必须为NULL。 6）Engine *eng OpenSingleUse(const char *startcmd,void *dcom,int *retstatus);

启动一个只允许用户使用的MATLAB引擎，在Windows系统中startcmd和dcom参数始终为NULL，retstatus为返回engOpenSingleUse函数的执行状态。

7）int engOutputBuffer(Engine *ep,char *p,int n);

设置MATLAB引擎的输出内存，存储engEvalString函数执行字符串后的输出结果。其中n表示设置的输出内存可以存放的字符个数，如果输出结果的字符串大于n，则只存储前n个字符。

8）int engPutVariable(Engine *ep,const char *name,const mxArray *mp); 向MATLAB引擎工作空间中写入一个MATLAB阵列变量，其中name为写入变量在MATLAB引擎工作空间中的名字。

9）int engSetVisible(Engine *ep,bool value);

设置MATLAB引擎工作窗口是否可见的属性，如果value=true,则MATLAB引擎窗口可见；如果value=false,则MATLAB引擎窗口不可见。

2、Visual C++6.0 MFC工程中调用MATLAB引擎

首先创建一个单文档的Visual C++6.0 工程matlabenginetest,其中view类的基类选为CFormView。

添加到CMatlabenginetestView中的界面元素及其ID如下表所列。

表4.1 matlabenginetest CmatlabenginetestView 控件列表

ID IDC_SrartEngine IDC_DrawSinc IDC_CloseEngine IDC_HIDEENGINECHECK IDC_CMDEDIT IDC_EVALUATESTRING IDC_OUTPUTEDIT 控件类型 Button Button Button 界面元素启动Engine 计算并绘制Sinc图象关闭Engine Check button 隐藏MATLAB引擎窗口 Edit Button Edit 命令输入编译框执行用户输入的命令 MATLAB引擎结果显示编译框 Matlabenginetest通过engEvalString函数调用MATLAB引擎执行相应的命令并计算显示sinc函数曲线，通过函数engOutputBuffer将m_outbuff设为MATLAB引擎输出的结果存储的缓冲区。另外需要注意的是，由于本工程调用MATLAB引擎的API函数，因此需要在工程设置中加入静态链接库libeng.lib。[7] 4.2.3 程序代码

/* * * * * * * * * * * * * * * * * matlabenginest.h * * * * * * * * * * * * * */

#if !defined(AFX_MATLABENGINETESTVIEW_H__2BA448E9_DCB7_4CAB_8DE3_46FED317C086__INCLUDED_) #define

AFX_MATLABENGINETESTVIEW_H__2BA448E9_DCB7_4CAB_8DE3_46FED317C086__INCLUDED_ #if _MSC_VER > 1000 #pragma once

#endif // _MSC_VER > 1000 #include \"engine.h\"

#define _MAX_BUFF_CHAR_NUM 2000 class CMatlabenginetestView : public CFormView {

protected: // create from serialization only

CMatlabenginetestView();

DECLARE_DYNCREATE(CMatlabenginetestView)

public:

//{{AFX_DATA(CMatlabenginetestView)

enum { IDD = IDD_MATLABENGINETEST_FORM }; CEdit m_outputedit; BOOL m_nIsCheck; CString m_cmdedit; //}}AFX_DATA

// Attributes public:

CMatlabenginetestDoc* GetDocument();

// Operations public:

Engine * m_ep;

char m_outbuff[_MAX_BUFF_CHAR_NUM];

// Overrides

// ClassWizard generated virtual function overrides //{{AFX_VIRTUAL(CMatlabenginetestView) public:

virtual BOOL PreCreateWindow(CREATESTRUCT& cs); protected:

virtual void DoDataExchange(CDataExchange* pDX); // DDX/DDV support virtual void OnInitialUpdate(); // called first time after construct

virtual BOOL OnPreparePrinting(CPrintInfo* pInfo); virtual void OnBeginPrinting(CDC* pDC, CPrintInfo* pInfo); virtual void OnEndPrinting(CDC* pDC, CPrintInfo* pInfo); virtual void OnPrint(CDC* pDC, CPrintInfo* pInfo); //}}AFX_VIRTUAL

// Implementation public:

virtual ~CMatlabenginetestView();

#ifdef _DEBUG

virtual void AssertValid() const;

virtual void Dump(CDumpContext& dc) const;

#endif protected:

// Generated message map functions protected: };

#ifndef _DEBUG // debug version in matlabenginetestView.cpp inline CMatlabenginetestDoc* CMatlabenginetestView::GetDocument() { return (CMatlabenginetestDoc*)m_pDocument; }

//{{AFX_MSG(CMatlabenginetestView) afx_msg void OnCloseEngine(); afx_msg void OnDrawSinc(); afx_msg void OnStartEngine(); afx_msg void OnPaint();

afx_msg void OnHideenginecheck(); afx_msg void OnEvaluatestring();

afx_msg void OnSize(UINT nType, int cx, int cy); afx_msg void OnChangeOutputedit(); //}}AFX_MSG

DECLARE_MESSAGE_MAP()

#endif

//{{AFX_INSERT_LOCATION}}

// Microsoft Visual C++ will insert additional declarations immediately before the previous line. #endif

// !defined(AFX_MATLABENGINETESTVIEW_H__2BA448E9_DCB7_4CAB_8DE3_46FED317C086__INCLUDED_)

/* * * * * * * * * * * * * * * matlabenginetest.cpp * * * * * * * * * * * * * * * * */ #include \"stdafx.h\"

#include \"matlabenginetest.h\" #include \"matlabenginetestDoc.h\" #include \"matlabenginetestView.h\" #ifdef _DEBUG

#define new DEBUG_NEW #undef THIS_FILE

static char THIS_FILE[] = __FILE__; #endif

// CMatlabenginetestView

IMPLEMENT_DYNCREATE(CMatlabenginetestView, CFormView) BEGIN_MESSAGE_MAP(CMatlabenginetestView, CFormView)

//{{AFX_MSG_MAP(CMatlabenginetestView) ON_BN_CLICKED(IDC_CloseEngine, OnCloseEngine) ON_BN_CLICKED(IDC_DrawSinc, OnDrawSinc) ON_BN_CLICKED(IDC_StartEngine, OnStartEngine) ON_WM_PAINT()

ON_BN_CLICKED(IDC_HIDEENGINECHECK, OnHideenginecheck) ON_BN_CLICKED(IDC_EVALUATESTRING, OnEvaluatestring) ON_WM_SIZE()

ON_EN_CHANGE(IDC_OUTPUTEDIT, OnChangeOutputedit)

//}}AFX_MSG_MAP // Standard printing commands

ON_COMMAND(ID_FILE_PRINT, CFormView::OnFilePrint)

ON_COMMAND(ID_FILE_PRINT_DIRECT, CFormView::OnFilePrint)

ON_COMMAND(ID_FILE_PRINT_PREVIEW, CFormView::OnFilePrintPreview)

END_MESSAGE_MAP()

// CMatlabenginetestView construction/destruction CMatlabenginetestView::CMatlabenginetestView() { }

CMatlabenginetestView::~CMatlabenginetestView() { }

void CMatlabenginetestView::DoDataExchange(CDataExchange* pDX) {

: CFormView(CMatlabenginetestView::IDD)

//{{AFX_DATA_INIT(CMatlabenginetestView) m_nIsCheck = FALSE; m_cmdedit = _T(\"\"); //}}AFX_DATA_INIT

// TODO: add construction code here m_ep = NULL;

memset(m_outbuff,0,_MAX_BUFF_CHAR_NUM*sizeof(char));

if(m_ep!=NULL) { }

engClose(m_ep);

CFormView::DoDataExchange(pDX);

//{{AFX_DATA_MAP(CMatlabenginetestView) DDX_Control(pDX, IDC_OUTPUTEDIT, m_outputedit);

}

DDX_Check(pDX, IDC_HIDEENGINECHECK, m_nIsCheck); DDX_Text(pDX, IDC_CMDEDIT, m_cmdedit); //}}AFX_DATA_MAP

BOOL CMatlabenginetestView::PreCreateWindow(CREATESTRUCT& cs) { }

void CMatlabenginetestView::OnInitialUpdate() { }

// CMatlabenginetestView printing

BOOL CMatlabenginetestView::OnPreparePrinting(CPrintInfo* pInfo) { }

void CMatlabenginetestView::OnBeginPrinting(CDC* /*pDC*/, CPrintInfo* /*pInfo*/) { }

void CMatlabenginetestView::OnEndPrinting(CDC* /*pDC*/, CPrintInfo* /*pInfo*/) {

// TODO: Modify the Window class or styles here by modifying // the CREATESTRUCT cs

return CFormView::PreCreateWindow(cs);

CFormView::OnInitialUpdate(); GetParentFrame()->RecalcLayout(); ResizeParentToFit();

m_outputedit.SetWindowText(\"\");

// default preparation

return DoPreparePrinting(pInfo);

// TODO: add extra initialization before printing

// TODO: add cleanup after printing

}

void CMatlabenginetestView::OnPrint(CDC* pDC, CPrintInfo* /*pInfo*/) { }

// CMatlabenginetestView diagnostics #ifdef _DEBUG

void CMatlabenginetestView::AssertValid() const { }

void CMatlabenginetestView::Dump(CDumpContext& dc) const { }

CMatlabenginetestDoc* CMatlabenginetestView::GetDocument() // non-debug version is inline { }

#endif //_DEBUG

// CMatlabenginetestView message handlers void CMatlabenginetestView::OnCloseEngine() { }

void CMatlabenginetestView::OnDrawSinc()

// TODO: add customized printing code here

CFormView::AssertValid();

CFormView::Dump(dc);

ASSERT(m_pDocument->IsKindOf(RUNTIME_CLASS(CMatlabenginetestDoc))); return (CMatlabenginetestDoc*)m_pDocument;

// TODO: Add your control notification handler code here engClose(m_ep); m_ep=NULL;

{

// TODO: Add your control notification handler code here if(m_ep!=NULL) {

engEvalString(m_ep,\"x1=0.01:0.01:10*pi\"); engEvalString(m_ep,\"x2=-10*pi:0.01:-0.01\"); engEvalString(m_ep,\"y1=sin(x1)./x1\"); engEvalString(m_ep,\"y2=sin(x2)./x2\"); engEvalString(m_ep,\"plot([x2 x1],[y2 y1])\"); //engEvalString()中为调用的MATLAB程序

//产生一个错误，从而验证MATLAB引擎的输出已经被捕获

}

void CMatlabenginetestView::OnStartEngine() { }

void CMatlabenginetestView::OnPaint() {

}

engEvalString(m_ep,\"1/0\");

m_outputedit.SetWindowText(m_outbuff);

else { }

AfxMessageBox(\"请启动MATLAB引擎\

// TODO: Add your control notification handler code here m_ep = engOpen(NULL);

engOutputBuffer(m_ep,m_outbuff,_MAX_BUFF_CHAR_NUM);

CPaintDC dc(this); // device context for painting

// TODO: Add your message handler code here

}

// Do not call CFormView::OnPaint() for painting messages

void CMatlabenginetestView::OnHideenginecheck() { }

void CMatlabenginetestView::OnEvaluatestring() { }

void CMatlabenginetestView::OnSize(UINT nType, int cx, int cy) {

// TODO: Add your control notification handler code here if(m_ep==NULL) return; UpdateData(true);

engSetVisible(m_ep,!m_nIsCheck);

// TODO: Add your control notification handler code here UpdateData(true);

LPSTR pstr = m_cmdedit.GetBuffer(m_cmdedit.GetLength()); if(m_ep!=NULL) { } else { }

AfxMessageBox(\"请启动MATLAB引擎\engEvalString(m_ep,(const char*)pstr); m_outputedit.SetWindowText(m_outbuff);

CFormView::OnSize(nType, cx, cy);

// TODO: Add your message handler code here if(::IsWindow(m_outputedit.GetSafeHwnd()))

}

{ }

m_outputedit.MoveWindow(0,140,cx,cy-140,true);

void CMatlabenginetestView::OnChangeOutputedit() { }

// TODO: If this is a RICHEDIT control, the control will not

// send this notification unless you override the CFormView::OnInitDialog() // function and call CRichEditCtrl().SetEventMask() // with the ENM_CHANGE flag ORed into the mask. // TODO: Add your control notification handler code here

图4.2 VC调用MATLAB引擎的用户界面

10.80.60.40.20-0.2-0.4-40-30-20-10010203040

图4.3 通过调用方法实现绘制sinc图形

4．3 数据分组

本课题为一路数据采集，得到一组电压信号数据。

谱线的精确辨识需要相继两次DFT分析，考虑到完全对称窗的选取，需从数据采集得到的数据中取N+2个数据，这里取4098个数据{a1，a2，…, a4098 }。将a1和a4097的平均值作为第一个点，其余不变，得到第一组数据，记为{x（n）}，将a2和a4098的平均值作为第一点，其余依次，得到第二组数据，记为{y（n）}。两组数据相差一个采样周期。

4．4 加窗

4.4.1 窗函数的选取

对时域信号作DFT分析首先应作加窗处理，选取理想的窗函数的基本原则是：（1）窗函数频谱的主瓣应尽量窄,能量尽可能集中在主瓣内,从而在谱分析时获得较高的频率分辨力,在数字滤波器设计中获得较小的过渡带。

（2）窗谱的旁瓣高度应尽量小而且随频率尽快衰减,以减小谱估计时的泄漏失真,在设计数字滤波器时减小通带的波动,提高阻带的衰减。但主瓣既窄,旁瓣又小衰减又快的窗函数是不容易找到的,比如矩形窗的旁瓣很大,但其主瓣宽度是最窄的,因此,在数据处理时通常需要做综合考虑取其折中。

在应用窗函数时,除了要考虑窗谱本身的特性外,还应当充分考虑被分析信号的特点以及具体的处理要求。在谱分析中,如果仅要求分析谱峰的频率而不计及幅值的

精度,则可选用主瓣较窄而便于分辨的矩形窗。如果要分析窄带信号且具有较强的干扰噪声时,则应选用旁瓣幅度较小的窗函数如Hamming 窗等。如果干扰距离信号较远,则可以采用旁瓣衰减速度较快的窗函数如Papoulis 窗等。对于随时间衰减的信号如脉冲响应信号等,可以应用指数窗等来提高分析的信噪比[12]。

下面是几种常见的窗函数（1）矩形窗[5] 时域： wnRNn 窗谱： WRejWReN1j2

Nsin2幅频函数： WR sin2波形：

图4.4 矩形窗及其频谱

（2）三角形窗时域：

2nN1wn22nN10nN12N1nN12

窗谱：

N1N1jsin42j2WBeeN1sin2NN1sin4j2e2Nsin222

幅频函数：

Nsin24 WBNsin2（3）汉宁窗时域： wn窗谱：

WHejWHe幅频函数：

N1j2212n1cosRNn 2N1

22WH0.5WR0.25WRWR

N1N1当N>>1，N-1近似N，所以此式能使三部分之和相互抵消，能量更集中在主瓣，但主瓣宽度却比矩形窗宽一倍，为8N。

波形：

图4.5 汉宁窗及其频谱

（4）海明窗时域：

2nwn0.0.46cosRNn

N1窗谱：

WHMejWHMeN1j2

幅频函数：

2WHM0.WR0.23WRN波形：

2WRN 

图4.6 海明窗及其频谱

下面是五种常见窗函数的基本参数比较 [2] ：

表4.2 五种窗函数基本参数

窗函数类型矩形窗三角窗汉宁窗海明窗布莱克曼窗旁瓣峰值n/dB －13 －25 －31 －41 －57 过渡带宽度B 4/N 8/N 8/N 8/N 12/N 阻带最小衰减s/dB －21 －25 －44 －53 －74 一般电号主要含有整数次谐波，因而常采用基于余弦窗的组合窗，这类窗只要选取观测时间是信号周期的整数倍，其频谱在各次整数倍谐波频率处幅值为零，因而谐波之间不发生相互泄漏。即使信号频率作小范围波动，泄漏误差也较小。余弦窗一般可以表达为（H为窗的项数-1）

1wnNahcosh0H2nh 公式4.1 N窗的项数越多，主瓣宽度越大，从而引起频谱分辨力的降低。但同时较多项数的窗函数能够产生较大的旁瓣衰减，有利于提高频谱计算精度，但组合窗的项数一般不

大于4，容易推得

1wahcWRhWRh 公式4.2 h02H把式公式1代入公式2并简化得到

2sinHHajNhwsine•jah 公式4.3 hhsinh0h02NsinNN公式4.3表明，欲使组合窗具有线形相位特性，应满足ah0,H1故常用的窗为

h0H两项汉宁（Hanning）窗、三项布莱克曼（Blackman）窗、四项布莱克曼哈里斯（Blackman-Harris）窗。

如果测量时间为信号周期整数倍（且大于等于4），则基波和各次谐波频率归一化后正好位于上述组合窗幅频特性零点，此时整数次谐波相互干扰误差为零。当信号频率波动时，各次谐波频率在窗零点附近也作一定的波动，不同窗的零点附近特性不同，故产生的泄漏误差也不相同。相同条件下，窗宽L和采样点数同时增大，频谱泄漏减小。四项布莱克曼-哈里斯窗大多数谐波衰减90dB以上，性能最优，布莱克曼窗次之 [15] 。

4.4.2 布莱克曼——哈里斯窗（blackman harris）的定义

布拉克曼窗也称二阶升余弦窗。布拉克曼利用更多的矩形窗频谱线性组合构成布拉克曼窗，其偶对称表示为

NN2wnamcosmn n,,1,0,1

22Nm0K1单边表示

m2wn1amcosmn n0,1,2N1

Nm0K1单边表示的布拉克曼窗幅度频谱函数为

W1m0K1mam222WmWm RRNNK1这个窗函数中系数的选择应满足以下约束条件am1.0假如布拉克曼窗有K

m0 .

个非零的系数am,则其振幅频谱将由（2K-1）个中心频率不同的矩形窗频谱线性组合而成。显然，要使窗函数频谱的主瓣宽度窄，则K值不能选得很大。

哈利斯曾找到了布拉克曼窗的一组三项及四项非零系数，如表4.3所示。

表4.3 am系数值

最大副瓣电平 A0 A1 A2 A3 三项系数 -67dB 0.42323 0.49755 0.07922 四项系数 -92dB 0.35875 0.48829 0.14128 0.01168 由这些系数所构成的窗函数成为布拉克曼——哈利斯窗，可获得很小的副瓣电平。布拉克曼——哈利斯窗取三项时得最大副瓣电平为-67dB，取四项系数时可得最大副瓣电平为-92dB [4] 。

在本课题设计中，我采用的是四项布莱克曼——哈里斯窗，其窗函数为：

其波形图如图4.7所示：

(a)

(b)

图4.7 布莱克曼-哈里斯窗及其频谱

a)时域特性 b)频域特性

4.4.3 窗函数的使用

设连续时域信号为xt，截断用窗函数为wt，其宽度为T，则相应的分辨力为△f=1/T，其傅立叶变换谱为Wf，即

Wf = Fw(t)

周期为Ts的采样脉冲序列gt可写为

gt =

ntnTs

其傅立叶变换谱为

Gf = fsm(fmfs)

式中fs=1Ts为采样频率, 则相应的数字窗wnTs 的窗谱WGf 为

WG(f) = Fw(t)g(t) = Wf * Gf

nN1N2w(nTs)ej2fnTs

= fsmWfmfs

对信号xt的DFT分析实质上是xt的加权截断采样脉冲序列

xwg(t) = xtwtgt 公式4.4

的傅立叶变换在诸离散频率点f=k△f上的谱密度采样值，即

X(k) = F(xwg(t))fkf= XWG(f)fkf 公式4.5

设信号xt具有有限条分立谱线,即

X(f) = Aieji(ffi) 公式4.6

i则依公式4.4可得

XWG(f) = X(f) * WG(f)

= AiejiWG(ffi)

i引用公式4.5、公式4.6有

X(k) = AiejiWG(kffi) 公式4.7

i虽然公式4.7表明X(k)是具有真实谱线的函数，但如果采用具有实域化窗谱的窗函数，以致WG(f)在主瓣之外可完全忽略。则当诸真实谱线频差大于WG(f)的主瓣宽度B，即fifj > B (ij) 时，只可能有一条与kf相差在主瓣半宽范围内

的真实谱线fi对X(k)有影响，即当kf-fi < B2时有X(k)  AiejiWG(kffi)

窗函数的频谱由一个主瓣和几个旁瓣组成，主瓣以时域信号的每个频率成分为中心。旁瓣在主瓣的两侧以一定的间隔衰减至零。FFT产生离散的频谱，出现在FFT每个谱线的是在每个谱线上的连续卷积频谱。如果原始信号的频谱成份与FFT中的谱线完全一致，这种情况下采样数据的长度为信号周期的整数倍，频谱中只有主瓣。没有出现旁瓣的原因是旁瓣正处在窗函数主瓣两侧采样频率间隔处的零分量点。如果时间序列的长度不是周期的整数倍，窗函数的连续频谱将偏离主瓣的中心，频率偏移量对应着信号频率和FFT频率分辨率的差异，这个偏移导致了频谱中出现旁瓣，所以，窗函数的旁瓣特性直接影响着各频谱分量向相邻频谱的泄漏宽度。

每种窗函数都有其自身的特性，选择正确的窗函数，必须先估计信号的频谱成分。如若信号中有许多远离被测频率的强干扰频率分量，应选择旁瓣衰减速度较快的窗函数；如果强干扰频率分量紧邻被测频率时，应选择旁瓣峰值较小的窗函数；如果被测信号含有两个或两个以上的频率成份，应选用主瓣很窄的窗函数；如果是单一频率信号，且要求幅度精度较高，则推荐用宽主瓣的窗函数。

虽然不存在副瓣的窗函数难以寻觅，但在本课题中采用布来克曼-哈利斯窗已能达到极高的精度，此数字窗为

w(n) = amcos((2N)mn) n = -N2，…，N2

m0k1其窗谱为

WG(f) = am[WR(fmf)WR(fmf)] 公式4.8 m02k1式中WR(f)为矩形窗谱

NN1n1wr(n) =  (N为偶数) 220其它其相应的窗谱为

WR(f) =

sin(NfTS)jfTSe

sin(fTS)将其代入公式4.8将使布莱克曼-哈里斯窗不具备实谱特性，因此有必要将数字矩窗修改为完全对称式

n1wr(n) = 0.5n0nN2N2 （N为偶数） N2此时

WR(f) =

sin(NfTS)

tg(fTS)这里布莱克曼-哈里斯窗选用四项系数窗（k=4），此时其最大副瓣电平相对于主瓣已衰减92dB，满足应用要求[1]。

在实际使用中，还可以使用高斯窗，高斯窗的定义如下：任何一种信号时宽T、频宽W之积TW满足下式

TW1 4在最小时宽频宽积（即TW1）的优化准则下得出的是高斯波形，因此，可4选择高斯波形作窗函数。高斯波形是一直延伸到无穷大去的，但是若在波形三倍均方值的地方进行截断，则误差很小。高斯窗如下式所示：

wne1n2N220nN1 2有限列长的高斯形窗相当于高斯形窗与矩形窗相乘。因此，有限列长的高斯形窗的频谱函数应是高斯窗频谱与矩形窗频谱的卷积，即

Wn12e2a122*WR2N2e212

当a2.5,很小时a是标准差值的倒数，是频谱宽度的一种度量[2]。

4．5 DFT变换

要了解快速傅立叶变换（简称FFT），首先必须对离散傅立叶变换有充分的理解。

knXkDFTxnxnWN k = 0，1，…，N1

n0N1因直接计算DFT的计算量与变换区间长度N的平方成正比，当N较大时，计算量太大，所以在快速傅立叶变换FFT出现之前，直接用DFT算法进行谱分析和信号的实时处理是不切实际的。FFT算法就是不断地把长序列的DFT分解成几个短序列的DFT，并利用WNkn的周期性和对称性来减少DFT的运算次数。FFT将DFT算法的运算效率提高了1到2个数量级。

FFT算法基本上分为两大类：时域抽取法FFT（DIT-FFT）和频域抽取法FFT（DIF-FFT）。

4.5.1 按时域抽取算法（DIT-FFT）

序列x（n）的N（N＝2M）点DFT为[2]

X(k)DFTxnN点xnWn＝0N－1knNk0,1,,N1

将上面的和式按n的奇偶性分解为

XkNn偶数2lxnWknNn奇数N2lxnWknN

k2l1Nk2l =x2lWNl0x2l1Wl0令x1lx2l，x2lx2l1。

2klkl因为WNWN2，所以上式可写成

N21l0N21XkxlW1klN2WkNxlW2l0klN2k0,1,,N1

按n的奇偶性将x（n）分解为两个N/2长的序列x1l和x2l，则N点DFT可分解为两个N/2点DFT来计算。用X1k和X2k分别表示x1l和x2l的N/2点DFT，即

X1kDFTx1lN/2点X2kDFTx2lN/2点N21l0x1lWNkl2k0,1,,N1 2N1 2N21l0klx2lWN2k0,1,,将式和式代入式，并利用WNkN2WNk和X1k、X2k的隐含周期性可得到：

kXkX1kWNX2kNkXkX1kWNX2k2k0,1,,N1 2这样，就将N点DFT的计算分解为计算两个N/2点离散傅立叶变换X1k和

X2k，再计算上式。为了将如下分解过程用运算流图表示，以便估计其运算量，观察运算规律，总结编程方法，先介绍蝶形图。

蝶形图及其运算功能如图4.8所示

A+CB

A-CB

图4.8 蝶形运算图

当N=23=8时，用蝶形图表示第一次时域奇偶抽取分解及上式的运算流图如图4.9所示

根据图4.9可以求得第一次分解后的运算量。图4.9包括两个N/2点DFT和N/2

NNNN个蝶形，每个点DFT需要次复数乘法和1次复数加法运算，每个蝶

2222形只有一次复数乘法运算和两次复数加法运算。所以，总的复数乘法次数为

2NNN2N 2N1N12222X(0) X(2) X(4) X(6) X(1) X(3) X(5) X(7) N/2点 DFT N/2点 DFT X1(0) X1(1) X1(2) X1(3) X2(0) X2(1) X2(2) X2(3) 0 WN1 WN2X(0) X(2) X(4) X(6) X(1) X(3) X(5) X(7) WN23 WN 图4.9 8点DFT一次时域抽取分解运算流图

总的复数加法次数为

NN2NN 1222222由此可见，经过一次抽取分解，当N>>1时，使N点DFT的运算量近似下降一半。所以，应当继续分解下去。将每个N/2点DFT分解为2个N/4点DFT，……，依次类推，经过M级时域奇偶抽取，可分解为N个1点DFT和M级蝶形运算，每级有N/2个蝶形。而1点DFT就是1点时域序列本身。N=8点DIT-FFT的运算流图如图4.10所示。

x(0) x(4) x(2) x(6) x(1) x(5) x(3) x(7) A(0) A(1) A(2) A(3) A(4) A(5) A(6) A(7) 0WN 0WN A(0) A(1) A(2) A(0) A(1) A(2) A(0) X(0) X(1) X(2) X(3) W A(3) 0NW A(4) A(5) 0NW 2NA(3) A(4) A(5) A(6) 0WN 1WN X(4) X(5) X(6) A(7) X(7) W A(6) A(7) 0WN 0NWN2 W 3NW2N A(7) 图4.10 8点DIT-FFT运算流图

4.5.2 按频域抽取算法（DIF-FFT）

频率推导的方式是将频域的元素重新作排列，找出与时域数据之间的关系。设序列x(n)长度为N=2M，首先将x(n)前后对半分开，得到两个子序列，其DFT可表示为如下形式：

kn X(k) = DFT[x(n)] = x(n)WNn0N1N21=

x(n)Wn0knN +

nN2N21knNx(n)Wx(nN1knN

N21 =

x(n)Wn0 +

n0Nk(nN2))WN 2N21 =

n0NknkN2x(n)Wx(n)WN 公式4.9 N2式中WNkN2 = (1)k，k为偶数时等于1，k为奇数时等于-1。

将X(k)分解成偶数组与奇数组，当k取偶数（k=2r，r=0，1，…，N21）时

N21= X(2r)

n0N2rnx(n)x(n)WN 2 .

N21 =

n0Nrnx(n)x(n)WN2 公式4.10 2当k取奇数（k=2r1，r=0，1，…，N21）时，

N21X(2r1) =

n0Nn(2r1)x(n)x(n)WN 2Nnnrx(n)x(n)WNWN2 公式4.11 2N21=

n0N)2令 ，n = 0，1，2，…，N21

Nng(n)x(n)x(n)WN2f(n)x(n)x(n将f(n)和g(n)分别代入4.10和4.11式，可得

N21rnX(2r)f(n)WN2n0 公式4.12 N21rnX(2r1)g(n)WN2n0公式4.10表明，X(k)按奇偶k值分为两组，其偶数组是f(n)的N2点的DFT，奇数组则是g(n)的N2点的DFT。

这种算法是对X(k)进行奇偶抽取分解的结果，所以称之为频域抽取法FFT。观察图4.11可知，DIF-FFT算法共M级运算，每级共有N2个蝶形运算，输入是序列为自然顺序，而输出为倒序排列。因此，M级运算完后，要对输出数据进行倒序才能得到自然顺序的X(k)。

在本课题设计中，我采用的是用C++调用MATLAB引擎的方法，在MATLAB中有直接的fft函数，使用十分方便。

x(0) x(1) x(2) x(3) x(4) x(5) x(6) x(7) 0WN 1WN X(0) 0WN 0WN X(4) X(2) WN2 0WN X(6) X(1) 0WN 0WN X(5) X(3) WN2 W 3NWN2 0WN X(7) 图4.11 N=8时DIF-FFT的运算流图

4．6 用查询方式找出真实谱线

4.6.1 谱线辨识

考虑到FFT算法通常以下式进行

X(k) = x(n)w(n)en0N1j2nkN

因此对应于双边对称算法

N21nN2x(n)w(n)ej2nkN

在找出真实谱线之前应将X(k)改造成

X'(k)=X(k)*(1)k

由于采用具有实谱特性的窗函数对信号进行DFT分析，所以在窗谱的主瓣内将含有若干个由窗泄漏造成的相位相同的谱分析值。只要数字窗谱WG(f)为实谱，则主瓣内含有i段相位（）相同的X(k)分析值。据此可以判定，这些相邻的X(k)值源于频率为fi的同一条真实谱线。故有i条真实谱线。本课题需找出十条真实谱线。

令k= 0到2048，先找出一些相位基本一致的点，然后找出其中幅度最大的点，然后跳过这些相位一致的点及旁瓣中的点，继续往下找相位一致的点，找出其中幅度

最大的点，共找10个这样的点。

因为有两组数据x（n）、y（n），所以有两组各10个这样的点。 4.6.2 频率和幅度

计算频率、幅度时，先分析信号x(t)的时移信号

y(t) = x(tTS)

由傅立叶变换的时移特性可以写出

Y(f) = X(f)ej2fTs Aieji(ffi)ej2fTs

可导出过程可得

Y(k) = Aiej(i2fiTs)WG(kffi)

此时，利用Y(k)与X(k)的幅角差值可以判定频率fi的准确值

i = ArgY(k)-ArgX(k) = 2fiTs

fi =

i1ifs =

2Ts2设fi = lif，li =

iN 2在求出准确频率lif之后，可以方便地确定该谱线的幅度

U(i)=

X(k)WG((kli)f)

至此，真实谱线的频率、幅度已可以全部求出。

4.7数据处理程序代码

%导入数据

a=importdata('data.txt');实际导入数据的方法 N=4096; %数据分组 x1=a;

x1(1)=(a(1)+a(N+1))./2;

for i=1:N y1(i)=a(i+1); end

y1(1)=(a(2)+a(N+2))./2; %加窗 for n=1:N

w(n)=0.35875+0.48829*cos(2*pi*n./N)+0.14128*cos(2*pi*2*n./N)+0.01168*cos(2*pi*3*n./N)

end for n=1:N

x1(n)=x1(n).*w(n); y1(n)=y1(n).*w(n); end

%做DFT变换 x2=fft(x1,N); y2=fft(y1,N); %改造 for k=1:N

x2(k)=x2(k).*(-1)^(k+1); y2(k)=y2(k).*(-1)^(k+1); end %相位 x3=angle(x2); y3=angle(y2);

%查询方式寻找真实谱线

%基频 fs=5000;

j=int16(50*N/fs+0.5); %取整，四舍五入 k(1)=0; for i=j-3:j+2

if norm(x2(i))>=norm(x2(i-1)) & abs(angle(x2(i))-angle(x2(i-1)))<=1e-3 %找出幅值最大，相位差误差小于0.001（则认为相位相同）的谱线， else k(1)=i; end end

if k(1)==0 A(1)=0; else ts=1./fs;

xwc=y3-x3; %相位差 fo=xwc./(2*pi)*fs; %频率 fd=fs./N; l=fo./fd; fu=k.*fd;

%计算相应的窗谱 fu11=fu-fd; fu12=fu+fd; fu21=fu-2*fd; fu22=fu+2*fd; fu31=fu-3*fd; fu32=fu+3*fd;

WR0=sin(N*pi*fu*ts)./tan(pi*fu*ts);

WR11=sin(N*pi*fu11*ts)./tan(pi*fu11*ts); WR12=sin(N*pi*fu12*ts)./tan(pi*fu12*ts); WR21=sin(N*pi*fu21*ts)./tan(pi*fu21*ts); WR22=sin(N*pi*fu22*ts)./tan(pi*fu22*ts); WR31=sin(N*pi*fu31*ts)./tan(pi*fu31*ts); WR32=sin(N*pi*fu32*ts)./tan(pi*fu32*ts); %布莱客曼-哈里斯窗窗谱

WG=0.35875*WR0+0.244145*(WR11+WR12)+0.070*(WR21+WR22)+0.00584*(WR31+WR32);

A(1)=norm(x2)./WG; %幅度 end

for i=2:10

j=int16(l(1)*i+0.5); k(i)=0; for m=j-3:j+2

if norm(x2(m))>=norm(x2(m-1)) & abs(angle(x2(m))-angle(x2(m-1)))<=1e-3 k(i)=m; end if k(i)==0 A(i)=0; else

fo(i)=xwc./(2*pi)*fs; fd=fs./N; l(i)=fo./fd; fu=k.*fd; fu11=fu-fd; fu12=fu+fd; fu21=fu-2*fd; fu22=fu+2*fd;

fu31=fu-3*fd; fu32=fu+3*fd;

WR0=sin(m*pi*fu*ts)./tan(pi*fu*ts); WR11=sin(m*pi*fu11*ts)./tan(pi*fu11*ts); WR12=sin(m*pi*fu12*ts)./tan(pi*fu12*ts); WR21=sin(m*pi*fu21*ts)./tan(pi*fu21*ts); WR22=sin(m*pi*fu22*ts)./tan(pi*fu22*ts); WR31=sin(m*pi*fu31*ts)./tan(pi*fu31*ts); WR32=sin(m*pi*fu32*ts)./tan(pi*fu32*ts);

WG=0.35875*WR0+0.244145*(WR11+WR12)+0.070*(WR21+WR22)+0.00584*(WR31+WR32);

A(i)=norm(x2)./WG end end end

结论

经过这次毕业设计，我对电网谐波有了一定的认识。现在电力系统中谐波的危害极大，已经引起人们的高度重视。所以供电部门、电力用户和设备制造商应合理规划电网。

本文是基于快速傅立叶的方法分析电网谐波。通过数据采集卡从电网中取得数据后，对信号及其时移连续两次FFT分析,从而可以得到谱线精确的幅度、频率值，进而可以求出电网的纹波系数。

参考文献

1 张处武.分立谱线的精确辨认方法.淮阴工业专科学校学报.1999(4):3-4 2 丁玉美，高西全.数字信号处理.西安：电子科技大学出版社，2005

3 罗安.电网谐波治理和无功补偿技术及装备.北京：中国电力出版社，2006 4 王世一.数字信号处理.北京：理工大学出版社，2004 5 钱同惠.数字信号处理.北京：机械工业出版社，2005 6 吴镇扬.数字信号处理.北京：高等教育出版社，2006 7 刘维.精通MATLAB与C/C++混合程序设计.北京：航空航天大学出版社，2005 8 李海涛，邓樱.MATLAB程序设计教程.北京：高等教育出版社，2005 9 王育坚.Visual C++面向对象编程教程.北京：清华大学出版社，2004 10 佟勇臣.C/C++程序设计基础.北京:科学出版社,2003

11 谭锋.VisualC++程序设计实训教程.北京:科学出版社,2003

12 郑星亮，程洁，魏任之.数字信号处理中的窗效应及窗函数的应用原则.1997(2):33-37

13 梁鸿生，王凯.VC6.0环境下调用MATLAB的方法研究.西安理工大学学报2004（20）：78-81

14 杨斌文.浅谈电网谐波污染.学习园地.2006（1）：18-19 15 李红伟，李在玉.FFT分析电力系统谐波的加窗插值算法.电工技术杂志.2004（10）：61-63

致谢

本文是在张处武教授的悉心指导和自始自终的热情鼓励下完成的，在三个多月的毕业设计过程中，老师兢兢业业的工作作风、平易近人的品格、严谨的治学态度和诲人不倦的学者风范，使我深受教益，值此论文完稿之际，谨向老师表示最忠心的感谢和深深的敬意。

同时感谢杨定礼老师在毕业设计过程中给予我的帮助与支持。

在毕业设计的过程中，我有幸成为充满学术气氛和友爱互助氛围的集体中的一员，从刘佳同学、冯俪元同学、李捷同学的讨论与交流中，我受益非浅。并且在查阅资料及毕业论文撰写过程中，吴鹏翔同学、陈辉同学给予了我极大的帮助。我从同学们的关心与帮助中所汲取的教益及和他们结下的友谊，是不能用简单的语言来包容的。在此一并表示感谢。

向所有关心、支持和帮助过我的朋友们表示最诚挚的谢意。

因篇幅问题不能全部显示，请点此查看更多更全内容

查看全文