 Разработка, производство и продажа радиоэлектронной аппаратуры  |
Требуется программист в Зеленограде
- обработка данных с датчиков; ColdFire; 40 тыс. e-mail:jobsmp@pochta.ru |
#define _N_SAMPLES_ (512) // 32,64,128,256,512,1024 points
#define _F_SAMPLES_ (100.0) // Частота тестового сигнала, Hz
#define _A_SAMPLES_ (100.0) // Амплитуда тестового сигнала
#if !defined M_PI
#define M_PI (3.14159265359)
#endif
float x[ _N_SAMPLES_ ], y[ _N_SAMPLES_ ], F, A, Re, Im, Fi;
//---------------------------------------------------------------
// I=1 --> БПФ
// I=-1 ---> Oбратное БПФ
void BPF( float *x, float *y, int N, signed char I )
{
volatile float c, s, t1, t2, t3, t4, u1, u2, u3;
volatile int i, j, p, l, L, M, M1, K;
L = N;
M = (N / 2);
M1 = (N - 1);
while( L >= 2 )
{
l = (L / 2);
u1 = 1.0000;
u2 = 0.000;
t1 = M_PI / (float)l;
c = cos( t1 );
s = (-1.000) * I * sin( t1 );
for( j = 0; j < l; j++ )
{
for( i = j; i < N; i += L )
{
p = i + l;
t1 = *(x+i) + *(x+p);
t2 = *(y+i) + *(y+p);
t3 = *(x+i) - *(x+p);
t4 = *(y+i) - *(y+p);
*(x+p) = (t3*u1) - (t4*u2);
*(y+p) = (t4*u1) + (t3*u2);
*(x+i) = t1;
*(y+i) = t2;
}/*for*/
u3 = (u1*c) - (u2*s);
u2 = (u2*c) + (u1*s);
u1 = u3;
}/*for*/
L /= 2;
}/*while*/
j = 0;
for( i = 0; i < M1; i++ )
{
if( i > j )
{
t1 = *(x+j);
t2 = *(y+j);
*(x+j) = *(x+i);
*(y+j) = *(y+i);
*(x+i) = t1;
*(y+i) = t2;
}/*if*/
K = M;
while(j >= K)
{
j -= K;
K /= 2;
}/*while*/
j += K;
}/*for*/
}/*BPF*/
//----------------------------------------------------
//----------------------------------------------------
//----------------------------------------------------
//----------------------------------------------------
//----------------------------------------------------
//----------------------------------------------------
/* Моделирование синусоидального тестового сигнала 1 */
/* Один период синуса за N измерений */
/* P - Массив выборок размерностью N */
/* A - Амплитуда (Частота _F_SAMPLES_ не имеет значения)*/
void Sinsignal_For_Test_One_Period( float *P, float A, int N )
{
volatile int i;
volatile float r, re, re1, im, im1;
for( i = 0; i < N; i++ )
{
*(P+i) = A * cos( 2.00 * M_PI *(float)i / (float)N );
}/*for*/
}/*Sinsignal_For_Test_One_Period*/
//----------------------------------------------------
//----------------------------------------------------
/* Моделирование синусоидального тестового сигнала 2 */
/* F периодов синуса за N измерений */
/* P - Массив выборок размерностью N */
/* F - Частота, Гц */
/* A - Амплитуда */
void Sinsignal_For_Test( float *P, float F, float A, int N )
{
volatile int i;
volatile float r, re, re1, im, im1;
re = cos((2.00 * M_PI * F ) / (float)N);
im = sin((2.00 * M_PI * F ) / (float)N);
re1 = A;
im1 = 0.00;
for( i = 0; i < N; i++ )
{
*(P+i) = re1;
r = re1;
re1 = (r * re) - (im1 * im);
im1 = (im1 * re) + (r * im);
}/*for*/
}/*Sinsignal_For_Test*/
//----------------------------------------------------
//----------------------------------------------------
__C_task main( void )
{
unsigned int j, N;
N = _N_SAMPLES_; // Points
F = _F_SAMPLES_; // Hz
A = _A_SAMPLES_; // Amplituda
////////////////////////////////
// Задаем входной сигнал
////////////////////////////////
// Пример 1. F периодов синуса за N выборок
//Sinsignal_For_Test( x, F, A, N ); // Генерация sin
// Пример 2. Один период синуса за N выборок
Sinsignal_For_Test_One_Period( x, A, N ); // Генерация sin
// Пример 3. Импульсный сигнал
// !!!!!! Сделать _N_SAMPLES_ = 32.
// _А_SAMPLES_ _F_SAMPLES_ не имеют значения
// for( j=0; j<8; j++ )
// {
// x[j] = 1.000;
// }/*for*/
// Остальные х[9...31] нули. Все y[0...31] нули
////////////////////////////////
// Вызов БПФ
////////////////////////////////
BPF( x, y, N, 1);
////////////////////////////////
////////////////////////////////
////////////////////////////////
for( j = 0; j < N/2; j++ )
{
/////////////// Расчет амплитуд
Re = *(x + j);
Im = *(y + j);
////////////////////////////////////////////////////////////////
A = sqrt( (Re * Re) + (Im * Im) );
if( j )
{
A = (A * 2.000);
}/*if*/
A = ( A / (float)N);
/////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////
/////////////// Расчет фаз
Fi = 0.000;
if( (A < (-0.001)) || (A > 0.001 ) )
{
if( (Re > (-0.000000001)) && (Re < 0.000000001 ) )
{
if( Im > 0.000000001 ) Fi = 1.5708;
if( Im < (-0.000000001) ) Fi = (-1.5708);
}/*if*/
else
{
Fi = atan2( Im , Re );
Fi = ((180.000 * Fi) / M_PI);
}/*else*/
}/*if*/
////////////////////////////////////////////////////
// Вывод на экран
///////////////////////////////////////////////////
printf(" A[%d] = %.3f, Fi = %.2f \n", j, A, Fi );
///////////////////////////////////////////////////
}/*for*/
}/*main*/
//##################################################
//##################################################
/////////////////////////
// Сплайн
/////////////////////////
Для нахождения значений функции в промежутках между узлами (точками измерений) удобно применять сплайн-интерполяцию, одним из достоинств которой является то, что ВЫЧИСЛЕННЫЕ по ней значения функции в узлах (точках измерения) гарантированно совпадут с ИЗМЕРЕННЫМИ значениями, что не всегда достижимо при "обычной" интерполяции полиномами (Ньютона,Лагранжа,Эверетта, где для повышения точности
применяются полиномы высокой степени, что само по себе дает погрешность, да и время
вычисления увеличивается. В случае сплайна достаточно полинома 3-ей степени).
Ниже приведена Си-функция для вычисления по сплайн-методу значения Y в промежутках между известными узлами (точками измерений).
Работать с ней очень просто (покажу на примере 6 узловых точек (MAXIMAL=6)):
Пусть в известные Вам 6 точек времени Вы намеряли 6 значений и хотите построить
на экране плавненький график.
1. Сохраните измеренные значения (вольты, омы, градусы, литры...) во float-виде в глобальном массиве Yy[1]...Yy[6].
(Обратите внимание, что здесь индекс массива идет не с нуля, а с 1. Это рудимент Бейсика, т.к. первоначально программа была написана на нем, и работала на Синклере и БК).
2. ВременА измерений сохраните во float-виде в глобальном массиве Xx[1]...Xx[6] (Конечно, каждое время Xx[i] должно соответствовать "своему" измерению Yy[i] ).
Пусть для определенности временА измерений 1, 2, 3, 4, 5, 6 мс. Промежутки между измерениями, вообще-то, могут быть неодинаковыми. Сплайну пох.
3. Ну, и запускайте функцию задавая ей в качестве аргумента интересующее Вас время.
"Просканируйте" таким образом Ваш временной диапазон - и выводите
график на экран прибора.
4. С чувством выполненного долга займите очередь в кассу предприятия за премией.
Примеры вызова:
Value1 = Spline( 1.5 ); // вычислит значение функции для 1.5 мс
........................
Value2 = Spline( 4.0 ); // даст то же значение, что Вы намеряли в узловой точке 4 мс
........................
Исходный код на Си:
////////////////////////////////////////////
#define MAXIMAL 6
/* Для изменения количества узлов изменять только значение MAXIMAL, но не трогая
констант 6.0 внутри функции */
/*******************************************/
float power3(float zx )
{
return (zx*zx*zx);
}/*power3*/
/*******************************************/
float Spline( float TimePoint )
{
unsigned char ii, jj, kk;
float dd, ee, hh, ff, pp, xxx, rrr, yyy;
float Ll[MAXIMAL+2], Rr[MAXIMAL+2], Ss[MAXIMAL+2], Mm[MAXIMAL+2];
dd = Xx[2] - Xx[1];
ee = ( Yy[ 2 ] - Yy[ 1 ] ) / dd;
for( kk=2; kk <= MAXIMAL-1; kk++ )
{
hh = dd;
dd = Xx[kk+1] - Xx[kk];
ff = ee;
ee = ( Yy[kk+1] - Yy[kk] ) / dd;
Ll[kk] = dd / ( dd + hh);
Rr[kk] = 1.0 - Ll[kk];
Ss[kk] = ( 6.0 * (ee - ff) ) / ( hh + dd );
}/*for*/
for( kk = 2; kk <= MAXIMAL-1; kk++ )
{
pp = 1.0 / ( Rr[kk] * Ll[kk-1] + 2 );
Ll[kk] = (-Ll[kk]) * pp;
Ss[kk] = ( Ss[kk] - Rr[kk] * Ss[kk-1] ) * pp;
}/*for*/
Mm[ MAXIMAL ] = 0.0;
Ll[ MAXIMAL-1 ] = Ss[ MAXIMAL-1 ];
Mm[ MAXIMAL-1 ] = Ll[ MAXIMAL-1 ];
for( kk = MAXIMAL-2; kk >= 1; kk-- )
{
Ll[kk] = Ll[kk] * Ll[kk+1] + Ss[kk];
Mm[kk] = Ll[kk];
}/*for*/
xxx = TimePoint;
ii = 0;
if( xxx > Xx[ MAXIMAL ] )
{
// экстраполяция вправо по оси времени
dd = Xx[MAXIMAL] - Xx[MAXIMAL-1];
yyy = dd * Mm[ MAXIMAL-1 ] / 6.0 + ( Yy[MAXIMAL] - Yy[MAXIMAL-1] ) / dd;
yyy = yyy * ( xxx - Xx[MAXIMAL] ) + Yy[MAXIMAL];
}/*if*/
else if( xxx <= Xx[ 1 ] )
{
// экстраполяция влево по оси времени
dd = ( Xx[2] - Xx[1] );
yyy = (((-dd) * Mm[2]) / 6.0) + (( Yy[2] - Yy[1] ) / dd);
yyy = (yyy * ( xxx - Xx[1] ) ) + Yy[1];
}/*else if*/
else
{
// интерполяция
do{ ii++; }while( xxx > Xx[ ii ] );
jj = ( ii - 1 );
dd = ( Xx[ii] - Xx[jj] );
hh = ( xxx - Xx[jj] );
rrr = ( Xx[ii] - xxx );
pp = dd * dd / 6.0;
yyy = ( ( Mm[jj]*(power3(rrr)) + Mm[ii]*(power3(hh)) ) / 6.0 ) / dd;
yyy = yyy+((( Yy[jj]-(Mm[jj]*pp))*rrr) + ((Yy[ii]-(Mm[ii]*pp))*hh) ) / dd;
}/*else*/
return yyy;
}/*Spline*/
/*********************************/
Легко заметить, что задавать значения времени можно и за пределами периода
измерений, т.е. можно проводить экстраполяцию. Но не ждите от такой
экстраполяции чудес, здесь это будет всего лишь прямая линия. Это линейная
экстраполяция или т.н. Задание Асимптотического Поведения.
//###################################################
