Поможем написать учебную работу
Если у вас возникли сложности с курсовой, контрольной, дипломной, рефератом, отчетом по практике, научно-исследовательской и любой другой работой - мы готовы помочь.

Предоплата всего

Подписываем
Если у вас возникли сложности с курсовой, контрольной, дипломной, рефератом, отчетом по практике, научно-исследовательской и любой другой работой - мы готовы помочь.
Предоплата всего
Подписываем
#include "stdafx.h"
#include <stdlib.h>
#include <windows.h>
#include <complex>
#include <time.h>
#include "mor.h"
Real Sg1[Max_TL];
Real Sg2[Max_TL];
Real X[Max_TL];
Real Power1[Max_TL];
Real Power2[Max_TL];
Real Covar [Max_TL];
ComplexType Scpl1[Max_TLC];
ComplexType Scpl2[Max_TLC];
ComplexType SinCos[Max_TLC];
ComplexType zmas[Max_TLC];
float powermax; int imax;
volatile unsigned long R = time(NULL)%10; // randomize
//----------------------------------------------------------------------------
float Rnd(float d)
{
return d*random(256)/256;
}
//----------------------------------------------------------------------------
void SetDirFFT(short DirectFFT)
{
DirFFT = DirectFFT;
}
//----------------------------------------------------------------------------
void Init_fourier(void)
{
SetDirFFT(TO_FREQ);
}
//----------------------------------------------------------------------------
void Expand_With_Zero(ComplexArr x, int k)
{
int i,n;
n = (2 << (k-1));
for(i=n+1; i<=2*n; i++)
{
(x+i)->Re = 0;
(x+i)->Im = 0;
}
}
//----------------------------------------------------------------------------
void CalcTrigonom(ComplexArr sincos, int k) // sincos д. меняться в процедуре и возвращаться
{
int n, m, i, j, l;
float arg, b;
n = (2 << (k-1));
l = i = 1;
while(l<n)
{
j = l << 1;
b = (float) pi/l;
for(m=1; m<=l; m++)
{
arg = (1-m)*b;
(sincos+i)->Re = cos(arg);
(sincos+i)->Im = sin(arg);
i += 1;
}
l = j;
}
}
//----------------------------------------------------------------------------
void FFT(ComplexArr x, ComplexArr sincos, int k)
{
int option, n, m, mm, nn, ic, i, j, l, ii;
float real_p2, imag_p2, real_p1, imag_p1;
option = (DirFFT == TO_FREQ) ? 1 : -1;
n = (2 << (k-1));
mm = 0;
nn = n - 1;
for(m=1; m<=nn; m++)
{
l = n;
while(mm+l>nn)
l >>= 1;
mm = l + (mm % l);
if (mm>m)
{
real_p2 = x[m+1].Re;
x[m+1].Re = x[mm+1].Re;
x[mm+1].Re = real_p2;
imag_p2 = x[m+1].Im;
x[m+1].Im = x[mm+1].Im;
x[mm+1].Im = imag_p2;
}
}
l = 1;
ic = 1;
while (l<n)
{
ii = l << 1;
for(m=1; m<=l; m++)
{
real_p1 = sincos[ic].Re;
imag_p1 = option*sincos[ic].Im;
i = m;
ic += 1;
while (i<=n)
{
j = i + l;
real_p2 = real_p1*x[j].Re - imag_p1*x[j].Im;
imag_p2 = real_p1*x[j].Im + imag_p1*x[j].Re;
x[j].Re = x[i].Re - real_p2;
x[j].Im = x[i].Im - imag_p2;
x[i].Re = x[i].Re + real_p2;
x[i].Im = x[i].Im + imag_p2;
i = i + ii;
}
}
l = ii;
}
}
//----------------------------------------------------------------------------
void Simultan_FFT(ComplexArr x, ComplexArr y, ComplexArr sincos, ComplexArr z, int k)
{
int i, n;
n = 2 << (k-1);
for(i=1; i<=n; i++)
{
(z+i)->Re = (x+i)->Re;
(z+i)->Im = (y+i)->Re;
}
FFT(z, sincos, k);
for(i=2; i<=n; i++)
{
x[i].Re = (z[i].Re + z[n+2-i].Re)/2;
x[i].Im = (z[i].Im - z[n+2-i].Im)/2;
y[i].Re = (z[i].Im + z[n+2-i].Im)/2;
y[i].Im = (z[i].Re - z[n+2-i].Re)/2;
}
x[1].Re = z[1].Re;
x[1].Im = 0;
y[1].Re = z[1].Im;
y[1].Im = 0;
}
//----------------------------------------------------------------------------
void Convolve(ComplexArr x, ComplexArr y, int k)
{
int i, sign, n;
float dummy;
n = 2 << (k-1);
sign = (DirFFT == TO_FREQ)? -1 : 1;
for (i=1; i<=n; i++)
{
dummy = x[i].Re*y[i].Re - sign*x[i].Im*y[i].Im;
x[i].Im = x[i].Im*y[i].Re + sign*x[i].Re*y[i].Im;
x[i].Re = dummy;
}
}
//----------------------------------------------------------------------------
void Power_spectrum(ComplexArr x, RealArr power, int k)
{
int i, n;
n = 2 << (k-1);
for (i=1; i<=n; i++)
power[i] = sqrt(x[i].Re*x[i].Re + x[i].Im*x[i].Im);
powermax=power[0];
for (i=1; i<=n; i++)
if (powermax<power[i])
{powermax=power[i];
imax=i;
}
}
//----------------------------------------------------------------------------
int random(int value)
{
#define Random(Max) ((R=(R*9301L+49267L)%233280L)%(long)Max) // random
return Random(value);
}
//----------------------------------------------------------------------------
void Model_EEG(RealArr sg1, RealArr sg2)
{
int ii,jj,k1;
int AEEG[3]={50,20,30}; //амплитуды для альфа, бетта и тетта-ритма
int FreqEEG[3]={10,20,5};//частоты для ритмов ЭЭГ
int FiEEG[3]={1,3,0};//сдвиг по фазе для ритмов ЭЭГ
Real Eps_AEEG[3], Eps_FreqEEG[3], Eps_FiEEG[3];
for (ii=0; ii<=2; ii++)
{
Eps_AEEG[ii]=AEEG[ii];
Eps_FiEEG[ii]=FiEEG[ii];
}
for (ii=1; ii<=Max_TL; ii++) // меняем условие пока ii меньше Max_TL делать ...
{
sg1[ii]=0;
for (jj=0;jj<=2;jj++)// поставить условие цикл работает пока jj<=2
{
k1=(AEEG[jj]+random(Eps_AEEG[jj]))*sin(2*pi*0.01*FreqEEG[jj]*ii+FiEEG[jj]+Eps_FiEEG[jj]);
sg1[ii]=sg1[ii]+k1;//моделируемые сигналы
}
}
for (ii=1; ii<=Max_TL; ii++)// меняем условие пока ii меньше Max_TL делать ...
{
sg2[ii]=0;
for (jj=0;jj<=2;jj++)// поставить условие цикл работает пока jj<=2
{
k1=+(AEEG[jj]+random(Eps_AEEG[jj]))*sin(2*pi*0.01*FreqEEG[jj]*ii+FiEEG[jj]+Eps_FiEEG[jj]);
sg2[ii]=sg2[ii]+k1; //моделируемые сигналы
}
}
}
//----------------------------------------------------------------------------
void Model_Rect(RealArr sg1, RealArr sg2)
{
int ii;
for (ii=1; ii<=Max_TL; ii++)
sg1[ii] = sg2[ii] = (ii<Max_TL/2)?50:-50;
}
//----------------------------------------------------------------------------
void Model_Sin(RealArr sg1, RealArr sg2)
{
int ii;
for (ii=1; ii<=Max_TL; ii++)
sg1[ii] = sg2[ii] = 100*sin(10*pi*ii/180);
}
//----------------------------------------------------------------------------
void Model_Saw(RealArr sg1, RealArr sg2)
{
int ii;
for (ii=1; ii<=Max_TL; ii++)
sg1[ii] = sg2[ii] = ii%100;
}
//----------------------------------------------------------------------------