|
|
Вот один из вариантов.
Взят из рабочей программы. Несколько упрощен.
Краткое объяснение.
Есть сигнал (signal), в котором кроме сигнала оказалась 50 Гц помеха.
Требовалось ее убрать. Достаточно даже было первую гармонику ее убрать. Но не искажать сигнал. Т.е. всякие переходные процессы исключались. Т.е. нужно было оценить частоту (ее уточнить), амплитуду и фазу. Затем сгенерировать cos с такими параметрами и вычесть его из сигнала. Генерация косинуса и вычитание выкинуты.
Использовалось окно Гегенбауэра 4 степени. Плюс дополнение 0 сигнального массива до учетверенной длины.
В программе lensign - длительность сигнальной реализации в секундах,
freqp - предполагаемая частота помехи.
Вариант приводимый ниже прост. Более сложный - аппроксимация спектра в требуемой области. Как амплитуды, так и фазы.
Вот только не знаю, что html оставит от программы.
color=7f7f7f>,4096);
for(i=0;i<4096;i++) {re[i]/=4096.;im[i]/=4096;}
i=(int)(floor((lensign*freqp)));
lsrh=(i-2)*4;
hsrh=(i+2)*4;
x=0.;
for(i=lsrh;i<=hsrh;i++)
{
y=re[i]*re[i]+im[i]*im[i];
if(y>x) {imax=i;x=y;}
}
sumsig=0.;
sumfr=0.;
for(i=(imax-2);i<=(imax+2);i++)
{
y=re[i]*re[i]+im[i]*im[i];
sumsig+=y;
sumfr+=y*(float)i;
}
amp=2.*sqrt((sumsig/sumgeg)); //amplitude estimation
frq=sumfr/sumsig; //frequency estimation
i=(int)(floor((frq+0.5)));
x=re[i];
y=im[i];
phase=atan2(y,x); //phase estimation
/* or:
i=(int)(floor((frq)));
x=atan2(im[i],re[i]);
y=atan2(im[i+1],re[i+1]);
phase=(y-x)*(frq-(float)i)+x; //phase estim.
*/
}
