基于FFT和LZW波形数据压缩.doc

《基于FFT和LZW波形数据压缩.doc》由会员上传分享，免费在线阅读，更多相关内容在教育资源-天天文库。

基于FFT的电能质量波形数据的LZW压缩算法第一章绪论1.1前言随着科技的进步，现代电力系统中用电负荷结构发生了重大变化。其非线性、冲击性以及不平衡等的负荷特性，使电网的电压波形发生畸变或引起电压波动和闪变以及三相不平衡，甚至引起系统频率波动等，对供电电能质量造成严重的干扰或污染。电网正面对越来越多的电能质量问题，这使得电能质量的研究十分紧迫。在电能质量监控方面，有两个趋势：其中之一是智能化，智能化旨在减轻人的劳动，能自动对电能质量问题进行识别和数据处理，从而实现全面的无人监控功能。另一个则是远程化。随着电力工业的发展和电网规模的扩大，供电部门和用户都迫切需要对较大量的监测点进行监控，然而各点的分散，距离远近不同，监测电能质量的问题也根据用户和电网的需要而各不相同。所以远程化就可以适应不同层次的监控要求，从而使电能质量的监控点能够分布到电网中的任何地方，并且具有良好的在线功能。而远程化必然带来的问题就是，监测点和监控站之间的通信问题以及大量的电能质量数据的传输问题都十分重要。因此对数据进行压缩是减少数据存储费用和提高性能的有效办法，具有重大的实际意义。电能质量分析仪表需要采集大量的三相电数据，从数据中解析出工业用电的各种参数，以便指导生产过程。而由于嵌入式设备本身的CPU处理能力、通信能力和存储能力所限，要求对采集到的大量数据先进行压缩然后再传输、存储和分析等，以缓解嵌入式设备的通信带宽低、存储容量小、计算能力弱等缺点带来的不便。1.2国内外研究近况 随着计算机技术的飞速发展，数据压缩作为海量信息存储和传输的支撑技术受到了人们的极大重视。目前，对于数据压缩算法的研究主要集中在数字图像处理、语音信号的分析等研究领域，也已经取得了显著的成果，而数据压缩技术在电力系统中的应用则相对较少。目前电能质量问题主要的分析方法可分为时域、频域和基于数学变换的分析方法三种。时域分析方法如差值法等，主要用于电能质量扰动信号的检测，方法快速简单，适合于检测电压凹陷、暂态脉冲等暂态电能质量问题，而不能检测如电压偏差、谐波、周期性陷波等稳态、周期性电能质量问题。频域分析方法主要用于电能质量中的谐波问题的分析，包括频谱分布、谐波潮流计算等。基于数学变换的分析方法主要有傅立叶变换方法、矢量变换方法以及近年出现的小波变换方法和人工神经网络分析方法等。傅立叶变换作为经典的信号分析方法具有正交、完备等许多优点，如FFT这样的快速算法。在频域，傅立叶变换具有较好的局部化能力，特别是对于那些频率成分比较简单的确定性信号，傅立叶变换很容易把信号表示成各频率成分的叠加和形式；但在时域中，傅立叶变换是对整个时间段的积分，没有局部化能力。由于小波分析方法具有良好的时-频局部化特性，它通过对不同的频率成分采用逐渐精细的采样步长，可以聚焦到信号的任意细节，能很好地处理微弱或突变信号，特别适合于对非稳态畸变波形问题进行分析。1.3本文任务本文针对提供的离线电能质量波形文件进行压缩，利用了有损压缩和无损压缩的结合；有损压缩部分是利用了对数据FFT之后的取阀值省去了一部分接近于0的数据，无损部分选择了对于重复性比较高的数据压缩效率高的LZW算法。压缩完后，再解压，计算谐波，谐波精度仍然要满足以下要求：表1压缩要求压缩率要求：保证精度的要求下压缩率越高越好。压缩速度：优先级较低。 第一章FFT2.1快速傅立叶变换傅立叶变换在频域具有良好的局部化能力，可以较好地描述信号的频率特性，适合用于谐波等电力系统稳态或稳态扰动现象的分析，如果信号为稳态扰动（如正弦波形、谐波等），采用快速傅立叶变换对信号进行压缩，得到相应的频率谱，记录了原始信号中不同频率成分的分量，于是这里选择了傅立叶变换没有选择对于暂态扰动有优势的小波变换。变换之后只要记住其中较大的系数，即可记录原始扰动数据中的大部分信息，达到压缩数据的目的，这里便成为有损压缩精度损失的来源。2.2离散傅立叶变换（DFT）DFT是信号分析与处理中的一种重要变换。因直接计算DFT的计算量与变换区间长度N的平方成正比，当N较大时，计算量太大，所以在快速傅里叶变换(简称FFT)出现以前，直接用DFT算法进行谱分析和信号的实时处理是不切实际的。直到发现了FFT算法，这种情况才发生了改变。当序列的点数不超过时，它的点定义为（1）反变换定义为（2）二者形式相似，快速算法的原理一样，这里先就其正变换进行讨论。令，当依次取为时，可表示为如下的方程组：（3）由上式可见，直接按照定义计算点序列的点时，每行含个复乘和个加，从而直接按定义计算点的总计算量为个复乘和个加。当 较大时，很大，计算量过大不仅耗时长，还会因字长有限而产生较大的误差，甚至造成计算结果不收敛。所谓快速傅里叶变换就是能大大减少计算量而完成全部点计算的算法。下面介绍两种经典的算法基于频域抽取和时域抽取的算法。2.2.1基于频域抽取基2FFT算法这里仅介绍基2算法，即是2的整次幂的情况。由定义（4）把分成两半，即和，代入（4）式得（5）由于（5）式两项又可合并为（6）当为偶数时，注意到，，（6）式变为（7）当为奇数时，，（6）式变为（8）这样就把一个点序列（）的点（）计算化成了两个点序列（和）的点（和）计算。由划分成 和的计算量为个加，即和和个乘，即由于算出的点，是的点（）中为偶数的那一半，由算出的则是为偶数的那一半，故需要把偶数的抽出来放在一起作为的（）输出，同时把奇数的抽出来放在一起作为的（）输出。由于是频域变量，故这种算法称为频域抽取的算法。接着，两个点仍可用上述方法各经个乘个加划分成两个点（同时还要做相应的频域抽取），从而共划分成4个点，总划分计算量仍是个加和个乘。这样的划分可一步步做下去，不难看出，每步的总划分计算量都是个加和个乘。经过步的划分后就划成了个如下两点的计算问题（9）上式所需计算量是2个加和1个乘，于是完成个两点的总计算量仍是个加和个乘。从而完成全部点的总计算量个加和个乘，这比直接按定义计算所需的个乘和加要少得多。例如，，，用算法计算所需的乘法个数为，而直接按定义计算所需的乘法个数为，二者相差倍。若直接计算需半小时，而用计算只需9s即可完成，可见其效率之高，而且越大，的效率提高越明显。 f[0]F[0]000F[0]F[0]000f[1]-1W20F[4]100F[2]F[1]001g[n]f[2]-1W40F[2]010F[4]F[2]010f[3]-1W41-1W20F[6]110F[6]F[3]011f[4]-1W80F[1]001F[1]F[4]100f[5]-1W81-1W20F[5]101F[3]F[5]101p[n]f[6]-1W82-1W40F[3]011F[5]F[6]110f[7]-1W83-1W41-1W20F[7]111F[7]F[7]111图1频域抽取的8点FFT计算流图一般情况下，由于做了次分奇偶的抽取，此算法最后的个两点计算出的不是顺序抽取的。次序的变化可用二进码来说明：第一次抽取所分的奇偶是由二进码第1位是1或0来区分的，该位为0时为偶数，该位为1时为奇数，第二次抽奇偶是由二进码第2位是1或0来区分的……，每次抽取都是把偶数项放在前（左）边，把奇数项放在后（右）边，从而抽取以后数的二进码是按照二进制位从左向右依次排列的，和普通二进制数从右向左依次排列的的规律正好相反，所以称为倒位序。在计算出之后要把倒位序变成顺序。所谓逆变换是指由求的计算，若直接按定义做计算，则除了求和号和正变换相同的计算量外，每算一个都还需再多做一个乘的乘法运算。故按定义完成全部点的总计算量是个加和个乘。下面从图导出它的快速算法，先讨论第3列的2点的逆运算如何完成。由式（7）得，由上式不难解出 （10）F[0]000F[0]F[0]0001/8f[0]F[1]001F[2]F[4]1001/8W2-0-1f[1]F[2]010F[4]F[2]0101/8W4-0-1f[2]F[3]011F[6]F[6]1101/8W2-0-1W4-1-1f[3]F[4]100F[1]F[1]0011/8W8-0-1f[4]F[5]101F[3]F[5]1011/8W2-0-1W8-1-1f[5]F[6]110F[5]F[3]0111/8W4-0-1W8-2-1f[6]F[7]111F[7]F[7]1111/8W2-0-1W4-1-1W8-3-1f[7]图2频域抽取的8点IFFT计算流图此计算过程如图2所示，可以看出：左边各列的划分计算也都是由个碟形运算来完成的，只是各碟形运算所乘的相移因子不同。把每个碟形运算都用图的办法变成对应的逆运算，并把它们按输入在左、输出在右重新排列，就得到了全部点的计算流图。给出了的示例，图中先对顺序输入的做次的频域抽取，并把个乘的运算合成了一个乘的运算放在了最前边，然后就开始做求逆的碟形运算。2.2.2基于时域抽取的基2FFT算法比较正变换和反变换的定义式可见，去掉乘的运算，把换成，交换、和、 ，反变换定义式就变成了正变换的定义式。对图2做这些变换，则得到图3的流程图。对图1做这些变换，则得到图4的流程图。这就是时域抽取的算法流图。进行碟形运算之前，先要对顺序的时域输入序列进行次的奇偶抽取，故称为时域抽取的算法。f[0]000f[0]f[0]000F[0]f[1]001f[2]f[4]100W20-1F[1]f[2]010f[4]f[2]010W40-1F[2]f[3]011f[6]f[6]110W20-1W41-1F[3]f[4]100f[1]f[1]001W80-1F[4]f[5]101f[3]f[5]101W20-1W81-1F[5]f[6]110f[5]f[3]011W40-1W82-1F[6]f[7]111f[7]f[7]111W20-1W41-1W83-1F[7]图3时域抽取的8点FFT计算流图这里先算出个两点的（11）f[0]1/8F[0]000F[0]F[0]000f[1]1/8-1W2-0F[4]100F[2]F[1]001f[2]1/8-1W4-0F[2]010F[4]F[2]010f[3]1/8-1W4-1-1W2-0F[6]110F[6]F[3]011f[4]1/8-1W8-0F[1]001F[1]F[4]100f[5]1/8-1W8-1-1W2-0F[5]101F[3]F[5]101f[6]1/8-1W8-2-1W4-0F[3]011F[5]F[6]110f[7]1/8-1W8-3-1W4-1-1W2-0F[7]111F[7]F[7]111 图4时域抽取的8点IFFT计算流图然后把个两点的组合成个4点的，再把个4点的组合成个8点的，经过次的组合之后，就得到了顺序点计算结果。图5就给出了用FFT算法与直接DFT算法所需乘法次数的比较曲线：图5FFT算法与直接DFT算法所需乘法次数的比较曲线2.3基2FFT算法编程以频域抽取的基2FFT正变换为例，对FFT的信号流图进行讨论，以找到FFT算法的规律。1）分级在进行变换的过程中，从点到两点共分了M级，如图1所示，从左到右依次为级，级，…，级。2）倒位序在频域抽取的基2FFT算法中，输出数据不是按照序列的先后顺序排列的，这是由于变换过程中，输出按奇、偶抽取的缘故。如果将序列中标号用二进制值表示，那么在FFT信号流图输入端，位于处，称为倒序。以8点FFT为例，顺序和倒序的关系如表2所示。 表2顺序和倒序对照表顺序倒序十进制数二进制数二进制数十进制数0123456700000101001110010111011100010001011000110101111104261537从表1可以看出，一个自然顺序二进制数，是在最低位加1，逢2向左移位；而倒序数的顺序是在最高位加1，逢2向右移位。用表示顺序数，表示倒序数，表示位权重。对于一个倒序数来说，下一个倒序数可以按下面的方法求：先对最高位加1，相当于十进制运算。如果，说明二进制最高位为0，则直接由得到下一个倒序值；如果，说明二进制最高位为1，则将的最高位变为0，通过实现，同时令，接着判断次高位是否为0，直到位为0时，令。归结起来算法流程图如图6所示。 输出j=N/2i=1toN-2t=f[i]f[j]=f[i]f[i]=f[j]k=N/2j=j-kk=k/2kjj=j+ki=k){j=j-k;k=k/2;}j=j+k;}{for(L=1;L<=m;L++)/*用蝶形运算完成FFT*/ {le=pow(2,L);B=le/2;/*同一蝶形结中参加运算的两个点的距离*/pi=3.14159;/*PI的取值也是进度损失的一部分*/for(j=0;j<=B-1;j++){p=pow(2,m-L)*j;ps=2*pi/N*p;w.real=cos(ps);/*W为系数商，当前系数与之前系数的商*/w.imag=-sin(ps);for(i=j;i<=N-1;i=i+le){ip=i+B;/*i,ip为参加运算的两个节点*/t=EE(xin[ip],w);xin[ip].real=xin[i].real-t.real;xin[ip].imag=xin[i].imag-t.imag;xin[i].real=xin[i].real+t.real;xin[i].imag=xin[i].imag+t.imag;}}}}return;}/*****************lzw算法***********************/#include#include#include#defineBITS12#defineHASHING_SHIFTBITS-8#defineMAX_VALUE(1<=next_code){*decode_stack=character;string=decode_string(decode_stack+1,old_code);}elsestring=decode_string(decode_stack,new_code);character=*string;while(string>=decode_stack)putc(*string--,output);if(next_code<=MAX_CODE){prefix_code[next_code]=old_code;append_character[next_code]=character;next_code++;}old_code=new_code;}printf(" ");/*cleanup...*/fclose(input);fclose(output);free(prefix_code);free(append_character);return1;}unsignedchar*decode_string(unsignedchar*buffer,unsignedintcode){ inti;i=0;while(code>255){*buffer++=append_character[code];code=prefix_code[code];if(i++>=4094){printf("Fatalerrorduringcodeexpansion. ");exit(0);}}*buffer=code;return(buffer);}unsignedintinput_code(FILE*input){unsignedintreturn_value;staticintinput_bit_count=0;staticunsignedlonginput_bit_buffer=0L;while(input_bit_count<=24){input_bit_buffer|=(unsignedlong)getc(input)<<(24-input_bit_count);input_bit_count+=8;}return_value=input_bit_buffer>>(32-BITS);input_bit_buffer<<=BITS;input_bit_count-=BITS;return(return_value);}voidoutput_code(FILE*output,unsignedintcode){staticintoutput_bit_count=0;staticunsignedlongoutput_bit_buffer=0L;output_bit_buffer|=(unsignedlong)code<<(32-BITS-output_bit_count);output_bit_count+=BITS;while(output_bit_count>=8){ putc(output_bit_buffer>>24,output);output_bit_buffer<<=8;output_bit_count-=8;}}/*****************mainprograme********************/#include#include#includeFILE*fp,*fp1,*fp2;typedefstruct_osc_eight{floatosc[8];/*osc[0]…osc[7]依次为同一时刻采样之Ua,Ub,Uc,Un,Ia,Ib,Ic,In*/}osc_eight;structcompxs[32696];intp[4096],q[4096];intNum=4096;constfloatpp=3.14159;main(intargc,char*argv[]){/*FFT部分*/osc_eighttmp_osc8;for(intj=0;j<8;j++)/*把八路数据依次导入进行FFT*/{fp=fopen("./ab.dat","rb");fp1=fopen("real.txt","a");/*建立不覆盖式real.txt文档*/fp2=fopen("imag.txt","a");/*建立不覆盖式imag.txt文档*/if(fp==NULL){perror("openab.dat");return-1;}inti=0;do{fread(&tmp_osc8,sizeof(osc_eight),1,fp);s[i].real=tmp_osc8.osc[j];s[i].imag=0;i++;}while(!feof(fp)); FFT(s,Num);for(intm=0;m<4096;m++){p[m]=(s[m].real/100);q[m]=(s[m].imag/100);fprintf(fp1,"%d,",p[m]*100);/*对数据取阀值后分别导出到建立的real,imag文档*/fprintf(fp2,"%d,",q[m]*100);}fclose(fp1);fclose(fp2);fclose(fp);}/*LZW部分*/printf("LZWcompressionanddecompressionutility ");if(argc!=4){printf(" Usage:lzw-c|dsourcefilenametargetfilename ");exit(0);}if(!strcmp(argv[1],"-c"))LZW_Compression(argv[2],argv[3]);/*进行压缩处理*/elseif(!strcmp(argv[1],"-d"))LZW_Decompression(argv[2],argv[3]);/*进行解压处理*/elseprintf(" Unknowcommand. ");return0;}