power.cpp

脑电信号分析软件

//---------------------------------------------------------------------------

#include <vcl.h>
#pragma hdrstop
#include "log.h"
#include "new.h"
#include "Power.h"
#include "Main.h"
#include "lvbo.h"
//---------------------------------------------------------------------------
#pragma package(smart_init)
#pragma resource "*.dfm"
#include <math.h>
#define FFTPOINT    512     //     256       1024
TForm_Power *Form_Power;
//---------------------------------------------------------------------------
void DrawRect(TCanvas *ptr,int left,int top, int right,int bottom, TColor color);
void Write_line(TCanvas *ptr, int x, int y, int x1 ,int y1, TColor color);
//---------------------------------------------------------------------------
 void DrawRect(TCanvas *ptr,int left,int top, int right,int bottom, TColor color)
{
    ptr->Brush->Color=color;
    ptr->FillRect(Rect(left,top,right,bottom));
}
//---------------------------------------------------------------------------
void Write_line(TCanvas *ptr,int x,int y, int x1,int y1, TColor color)
{
     if (x==x1 && y==y1)
     {
        ptr->Pixels[x][y]=color;
        return;
     }
	 ptr->Pen->Color=color;
	 ptr->MoveTo(x,y);
	 ptr->LineTo(x1,y1);
}
//---------------------------------------------------------------------------
void __fastcall TForm_Power:: DrawBackGroundGrid_PowerOut(int step,int q)//此处用于化功率输出的背景色及网格
{
    fs=100;
    int Offset_x=(Image2->Width-512)/2+10;
    int Offset_y=Image2->Height-220;
    int count1=q;
    int count2=0;
    int index_50= 512*2*50/fs;

    Image2->Canvas->Pen->Style=psDot;
    for(int i=Image2->Height-Offset_y;i>=0;i-=200/q)
        {
            Write_line(Image2->Canvas,Offset_x+20,50+i,Image2->Width-Offset_x+40,50+i,clSilver);//画横线
            Image2->Canvas->Brush->Style=bsClear;
            if (count1==0)
                Image2->Canvas->TextOut(Offset_x,50+i-5,"    0");      // 230
            else
                Image2->Canvas->TextOut(Offset_x,50+i-5,AnsiString(-count1*10));     // 230
            count1--;
        }
    Image2->Canvas->TextOut(Offset_x,50+0,"P/dB");     //   230

    for(int j=Offset_x;j<=512+Offset_x;j+=index_50)
        {
            Write_line(Image2->Canvas,j+20,50+0,j+20,50+Image2->Height-Offset_y,clSilver); //画竖线
            if (count2==0)
                Image2->Canvas->TextOut(j+20,50+Image2->Height-Offset_y+2,"0");
            else
                Image2->Canvas->TextOut(j-10+20,50+Image2->Height-Offset_y+2,AnsiString(count2*50));
            count2++;
        }
    Image2->Canvas->TextOut(512+Offset_x+25,50+Image2->Height-Offset_y,"f/Hz");
    Image2->Canvas->Pen->Style=psSolid;
}
//----------------------------------------------------------------------------
//----------------------------------------------------------------------------
 int  __fastcall TForm_Power:: Min_Max(double YComplex[])
{
    double q;
    double Min_Power=1;
    int quotient;
    int Max_i;
    int  NPoint =FFTPOINT;          //此处给出特定的值，

    Max_Power=0;
    for(int i=0;i<NPoint;i++)
    {
        q=YComplex[i];
        if (q<Min_Power)
            Min_Power=q;
        if (q>Max_Power)
        {
            Max_Power=q;
            //Max_i=i;
        }
    }
    //Label3->Caption=AnsiString(fs*Max_i/NPoint);
    quotient=-10*log10(Min_Power/Max_Power)/10+1;
    return quotient;
}
//----------------------------------------------------------------------------

void __fastcall TForm_Power:: Draw_PowerOut(int step,int q,double YComplex[])   //功率输出画图
{
    int i=0;
    int Offset_x=(Image2->Width-512)/2+30;
    for (int ix=0;ix<512 ;ix+=step)     //          Npoint/2
    {
        Write_line(Image2->Canvas,ix+Offset_x, 50+20-10*log10(YComplex[i]/Max_Power)*200/(q*10),ix+step+Offset_x,50+20-10*log10(YComplex[i+1]/Max_Power)*200/(q*10),clBlue);
        i++;
    }
}
//---------------------------------------------------------------------------

//---------------------------------------------------------------------------
__fastcall TForm_Power::TForm_Power(TComponent* Owner)
    : TForm(Owner)
{
}
//---------------------------------------------------------------------------
void FFTInit(double wx[],double wy[])   //初始化加权向量
{
	int i;
	double factr, angle;
	factr = 2.0 * M_PI / (double)FFTPOINT;

	for (i = 0;i < FFTPOINT;i++)
	{
		angle = (double)i * factr;
		wx[i] = cos(angle);
		wy[i] = sin(angle);
	}
}
int bitr(int k, int logn)
{
	short int ans = 0;
	short int j = k;

	for (short int i = 0;i < logn;i++)
	{
		ans = (ans << 1) + (j & 1);
		j = j >> 1;
	}
	return (ans);
}
void fft(double datax[],double datay[],double wx[],double wy[],int length)
{  // FFT运算函数，入口DATAX，DATAY，WX，WY  出口DATAX，DATAY
	short int n1, logn, i, j, k, l, logl, p;
	float sx, sy, tx, ty;   //   complex
	logn = 9;//10;
	n1 = length >> 1;
	j = logn - 1;
	k = 0;
	//------ transform

	for (logl = 0;logl < logn;logl++)
	{
		do
		{
			for (i = 0;i < n1;i++)
			{
				p = bitr((k >> j), logn);
				l = k + n1;
				//---------- conj
				sx = wx[p];
				sy = -wy[p];
				//---------- multiply
				tx = sx * datax[l] - sy * datay[l];
				ty = sx * datay[l] + sy * datax[l];
				//---------- sub
				datax[l] = datax[k] - tx;
				datay[l] = datay[k] - ty;
				//---------- add
				datax[k] = tx + datax[k];
				datay[k] = ty + datay[k];
				k++;
			} // dofor i

			k += n1;
            if(n1 == 0) break;   //////////////////////////////////////
		}
		while (k < length);

		k = 0;

		j--;

		n1 = n1 >> 1;
        if(n1 == 0) break;   /////////////////////////////////
	}
	//-------- reorder
	for (i = 1;i < length;i++)
	{
		k = bitr(i, logn);

		if (i > k)
		{
			//-------- exchange i,k elements
			sx = datax[i];
			sy = datay[i];
			datax[i] = datax[k];
			datay[i] = datay[k];
			datax[k] = sx;
			datay[k] = sy;
		}     // end if
	}   // end for
}
//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
void MagFFT(double Mag[],double datax[],double datay[],int length )
{    int i;
      double RealData,ImagData;
    for (i=0;i<length;i++)
    {
       RealData = (datax[i])*(datax[i]) ;
       ImagData=  (datay[i])*(datay[i]);
       Mag[i]= ((RealData+ImagData))/( length);    //
    }
}
//----------------------------------------------------------------------------
//----------------------------------------------------------------------------
/*void __fastcall TForm_Power::FormCreate(TObject *Sender)
{
     int p,q;
     int iFileHandle,i,Image_Heigh;
     int iBytesRead;
    // int p,q;
     double FFTMag[FFTPOINT+5];
     double max,min,swing;
     double wx[FFTPOINT+5],wy[FFTPOINT+5];   //计算FFT用
     double datax[FFTPOINT+5],datay[FFTPOINT+5];  // 计算FFT用
     iBytesRead=FFTPOINT;

    if (MainForm->End_x-MainForm->Start_x+1>512)
          MainForm->End_x=512+MainForm->Start_x;

    for(p=MainForm->Start_x;p<MainForm->End_x;p++)
    {
      for(q=0;q<16;q++)
        {
         Fileleng[p-MainForm->Start_x].td[q]=MainForm->buf[p*16+q] ;
        }
     }
    if (MainForm->End_x-MainForm->Start_x+1<512)
    {
       for (i=p-MainForm->Start_x+1;i<512;i++)
         {
         for(q=0;q<16;q++)
            {
             Fileleng[i].td[q]=MainForm->buf[p*16+q] ;
            }
         }
    }
     int s=2*512/FFTPOINT;//计算每点的间隔
     int Number;
     for(p=0;p<FFTPOINT;p++)
     {
       datax[p]=Fileleng[i].td[MainForm->Lead_Num];//ctg[p].td[8];
       datay[p] = 0;
     }
        Image1->Canvas->FillRect(Image1->ClientRect);
        for (i = 0;i < iBytesRead-1;i++)
        {
            Image1->Canvas->MoveTo(i+39, Image1->Height*0.6-datax[i]*0.03);  //      Image1->Height*0.8-
            Image1->Canvas->LineTo(i+1+39,Image1->Height*0.6- datax[i+1]*0.03); //   Image1->Height*0.8-
        }
 //-----进行FFT变换-------------------------------------------------------
        FFTInit(wx,wy);
        fft(datax,datay,wx,wy,iBytesRead);
        MagFFT(FFTMag,datax,datay,iBytesRead);
//------功率输出画图---------------------------------------------------------
        Number=Min_Max(FFTMag);
        DrawBackGroundGrid_PowerOut(s,Number);   //画出输出信号功率谱的背景
        Draw_PowerOut(s,Number,FFTMag);   //画出输出信号的功率谱
 } */
 //-----------------------------------------------------------------------------

//------------------------------------------------------------------------------
void __fastcall TForm_Power::Button2Click(TObject *Sender)
{
Form_Power->Close();
}
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------



     /*double wx[FFTPOINT+5],wy[FFTPOINT+5];   //计算FFT用
    double datax[FFTPOINT+5],datay[FFTPOINT+5];  // 计算FFT用
    int i;
    double factr, angle;
    int st,et;
	factr = 2.0 * M_PI / (double)FFTPOINT;
    for (i = 0;i < FFTPOINT;i++)
	{
		angle = (double)i * factr;
		datax[i] = cos(angle);
		datay[i] = 0;
	}
    for (i = 0;i < FFTPOINT-1;i++)
	{
        Image1->Canvas->MoveTo(i,datax[i]*50 + Image1->Height/2);
        Image1->Canvas->LineTo(i+1,datax[i+1]*50 + Image1->Height/2);
	}

    FFTInit(wx,wy);
    fft(datax,datay,wx,wy);

    for (i = 0;i < FFTPOINT-1;i++)
	{
        Image2->Canvas->MoveTo(i,datax[i]*10000000 + Image1->Height/2);
        Image2->Canvas->LineTo(i+1,datax[i+1]*10000000 + Image1->Height/2);
	} */

//}



//---------------------------------------------------------------------------

//    int iFileHandle,i,Image_Heigh;
//    int iBytesRead;
//    int p,q;
//    double FFTMag[FFTPOINT+5];
//    double max,min,swing;
//    double wx[FFTPOINT+5],wy[FFTPOINT+5];   //计算FFT用
//    double datax[FFTPOINT+5],datay[FFTPOINT+5];  // 计算FFT用
//    struct myCTG ctg[FFTPOINT+5];

//    int s=2*512/FFTPOINT;//计算每点的间隔
//    int Number;

/*    if(OpenDialog1->Execute())
    {
        iFileHandle = FileOpen(OpenDialog1->FileName, fmOpenRead);
        iBytesRead = FileRead(iFileHandle, ctg, FFTPOINT * sizeof(struct myCTG));
        FileClose(iFileHandle);
        iBytesRead = iBytesRead/sizeof(struct myCTG);
*/
//  for(p=0;p<iBytesRead;p++)
//    {

//             datax[p]=ctg[p].td[8];
//             datay[p] = 0;

//    }
//        Image1->Canvas->FillRect(Image1->ClientRect);
//        for (i = 0;i < iBytesRead-1;i++)
//        {
//            Image1->Canvas->MoveTo(i+39, Image1->Height*0.6-datax[i]*0.03);  //      Image1->Height*0.8-
//            Image1->Canvas->LineTo(i+1+39,Image1->Height*0.6- datax[i+1]*0.03); //   Image1->Height*0.8-
//        }
 //-----进行FFT变换-------------------------------------------------------
//        FFTInit(wx,wy);
//        fft(datax,datay,wx,wy,iBytesRead);
//        MagFFT(FFTMag,datax,datay,iBytesRead);
//------功率输出画图---------------------------------------------------------
//        Number=Min_Max(FFTMag);
//        DrawBackGroundGrid_PowerOut(s,Number);   //画出输出信号功率谱的背景
//        Draw_PowerOut(s,Number,FFTMag);   //画出输出信号的功率谱
//----------------------------------------------------------------

//    }
//}

//---------------------------------------------------------------------------
    /*    Image2->Canvas->FillRect(Image2->ClientRect);
        max = FFTMag[300];
        min = FFTMag[300];
        for (i = 400;i < 500;i++)
        {
            if(FFTMag[i] > max)
                max = FFTMag[i];
            if(FFTMag[i] < min)
                min = FFTMag[i];
        }
        swing =( max - min);

        Image_Heigh=Image2->Height ;
        Image2->Width =iBytesRead;     //    3 *
        Image2->Picture->Bitmap->Width = iBytesRead;  // 3 *    *3       *3
        for (i =0;i <= iBytesRead/2-1;i++)   //(i = 0;i < iBytesRead-1;i++)
        {
           //Image2->Canvas->MoveTo(i,(FFTMag[i] - min)* (Image2->Height/2/swing)+Image2->Height/100);
           // Image2->Canvas->LineTo(i+3,(FFTMag[i+1]- min)* (Image2->Height/2/swing)+Image2->Height/100);
             Image2->Canvas->MoveTo(i ,0.8*Image_Heigh-(FFTMag[i] - min)* (Image2->Height/2/swing)+Image2->Height/10);
            Image2->Canvas->LineTo(i+1 ,0.8*Image_Heigh-(FFTMag[i+1] - min)* (Image2->Height/2/swing)+Image2->Height/10);
        } */






//---------------------------------------------------------------------------