Міністерство освіти та науки України
Національний університет «Львівська політехніка»
Розрахункова робота
з курсу «Комп’ютерна обробка вимірювальної інформації»
на тему
«Фільтрація шумових сигналів із завадою синусоїдної форми»
Варіант №15
Завдання:
1. Задано 2 із 3 параметрів:
Anois — амплітуда шумового вимірювального сигналу;
Azav — амплітуда завади;
К — коефіцієнт пропорційності,
пов’язані співвідношенням
2. Задано розрядність АЦП, яким вимірюємо дискретні значення вимірюваного шумового сигналу, спотвореного завадою
NADC = 8 ... 24 розряди
3. Оскільки вимірюємо двополярну напругу, то визначаємо кількість дискрет за рівнем сигналу:
,
так, наприклад для 16-розрядного АЦП кількість дискрет
хоча повна кількість дискрет для цілого сигналу
.
4. Генеруємо вимірювальний шумовий сигнал з амплітудою Anois
DataNois[i] = Anois * (random(Nd)/(float) Nd – 0.5)*2.;
За заданою формулою згенеруються числа від –32768 до +32768 (якщо задати Nd = 65536), масштабуються до –0,5 ... +0,5, потім до –1 ... +1 і на останок будемо мати випадкові числа від –Anois до Anois.
5. Виводимо графік вимірювального шумового сигналу.
6. Генеруємо заваду заданої форми, заданої частоти (f), заданого зсуву фаз (() та кількості дискрет (n).
Вибираємо час вимірювання 1 секунда. Щоб на графіку відображалась кількість періодів пропорційна f. Так для f = 3 буде 3 періоди.
Крім того треба забезпечити дискретність за частотою для перетворення Фур’є з кроком 1 МГц.
krok = 2*M_PI*f/n;
for (i = 0 ...)
DataZav[i] = Azav*sin(i*krok+()
7. Вивести графік завади.
8. Додаємо дискретні значення вимірювального шумового сигналу та завади
Data[i] = DataNois[i] + DataZav[i];
9. Виводимо графік сумарного сигналу.
10. Робимо перетворення Фур’є для сумарного сигналу і виводимо графік спектру.
11. Визначаємо середнє значення гармонік спектру сигналу Gser (додаємо всі гармоніки та ділимо на їх кількість)
12. Знаходимо середньоквадратичне відхилення ( випадкової величини Gser.
13. Виходячи із закону Релея розподілу випадкової величини
за формулою статистичної функції розподілу
знаходимо таке значення G( для якого отримана ймовірність задовольняє умові
p(G) — задане.
14. Маючи G( і Gser можемо знайти коефіцієнт обмеження спектру для методу фільтрації за рівнем
15. Проводимо фільтрацію за рівнем, тобто присвоюємо усім гармонікам, які є більшими за рівень ( ( G( середні значення гармоніки спектру Gser.
16. Графічно відображаємо спектр після фільтрації
17. Робимо зворотне перетворення Фур’є та отримуємо вимірюваний шумовий сигнал без завад
18. Відображаємо графік вимірюваного шумового сигналу за часом.
Дані:
№
Anois
Azav
К
NADC
Форма завади
f, Гц
(
n
p(G)
14
—
2,0
0,7
18
cинус
8
30
1150
0,95
Код програми:
#include<graphics.h>
#include<iostream.h>
#include<complex.h>
#include<stdlib.h>
#include<stdio.h>
#include<conio.h>
#include<math.h>
#define n 500
#define pG 0.95
#define Nadc 14
#define fi 30
#define K 0.7
#define f 8
#define Azav 2
int main (void)
{
int gdriver = DETECT, gmode, errorcode;
initgraph(&gdriver, &gmode, "");
errorcode = graphresult();
if (errorcode != grOk)
{
printf("Graphics error: %s\n", grapherrormsg(errorcode));
printf("Press any key to halt:");
getch();
exit(1);
}
randomize();
long int i, j;
double Y[n], X[n], DataNois[n], DataZav[n], Data[n], Xg, Yg;
double Anois, h;
double Xp, Xk, Xgp=20, Xgk=620;
double Ygp=20, Ygk=460, Yp, Yk;
double a[n], b[n], AM, BM;
double Pgp=240, Pgk=20, Pk=0, Pp=0;
double ag[n], bg[n], G[n];
double Amax=0, Amin=0, Bmax=0, Bmin=0;
long double Nd;
Anois=K*Azav;
Nd=pow(2,Nadc-1);
h=2*M_PI*f/n;
Xp=0; Yp=0; Xk=0; Yk=0;
line(Xgp,Ygk,Xgk,Ygk);
line(Xgp,Ygp,Xgp,Ygk);
outtextxy(40,20,"DataNois");
for (i=0;i<n;i++)
{
X[i]=0+i*h;
Y[i]=Anois*(random(Nd)/(float)Nd-0.5)*2.;
DataNois[i]=Y[i];
if (X[i]>Xk) Xk=X[i];
if (X[i]<Xp) Xp=X[i];
if (Y[i]>Yk) Yk=Y[i];
if (Y[i]<Yp) Yp=Y[i];
}
for (i=0;i<n;i++)
{
Xg=Xgp+(X[i]-Xp)*(Xgk-Xgp)/(Xk-Xp);
Yg=Ygk-(Y[i]-Yp)*(Ygk-Ygp)/(Yk-Yp);
putpixel (Xg,Yg,2);
}
getch();
cleardevice();
Xp=0; Yp=0; Xk=0; Yk=0;
outtextxy(40,20,"DataZav");
line(Xgp,Ygk,Xgk,Ygk);
line(Xgp,Ygp,Xgp,Ygk);
for (i=0;i<n;i++)
{
X[i]=0+i*h;
Y[i]=Azav*sin(i*h+fi);
DataZav[i]=Y[i];
if (X[i]>Xk) Xk=X[i];
if (X[i]<Xp) Xp=X[i];
if (Y[i]>Yk) Yk=Y[i];
if (Y[i]<Yp) Yp=Y[i];
}
for (i=0;i<n;i++)
{
Xg=Xgp+(X[i]-Xp)*(Xgk-Xgp)/(Xk-Xp);
Yg=Ygk-(Y[i]-Yp)*(Ygk-Ygp)/(Yk-Yp);
putpixel (Xg,Yg,2);
}
getch();
cleardevice();
Xp=0; Yp=0; Xk=0; Yk=0;
outtextxy(40,20,"DataNois+DataZav");
line(Xgp,Ygk,Xgk,Ygk);
line(Xgp,Ygp,Xgp,Ygk);
for (i=0;i<n;i++)
{
X[i]=0+i*h;
Y[i]=DataNois[i]+DataZav[i];
Data[i]=Y[i];
if (X[i]>Xk) Xk=X[i];
if (X[i]<Xp) Xp=X[i];
if (Y[i]>Yk) Yk=Y[i];
if (Y[i]<Yp) Yp=Y[i];
}
for (i=0;i<n;i++)
{
Xg=Xgp+(X[i]-Xp)*(Xgk-Xgp)/(Xk-Xp);
Yg=Ygk-(Y[i]-Yp)*(Ygk-Ygp)/(Yk-Yp);
putpixel (Xg,Yg,2);
}
getch();
cleardevice();
for (j=0;j<n;j++)
{
a[j]=0; b[j]=0;
for (i=0;i<n;i++)
{
a[j]=a[j]+Data[i]*cos(2*M_PI*j*i/n);
b[j]=b[j]+Data[i]*sin(2*M_PI*j*i/n);
}
a[j]=a[j]/n;
b[j]=b[j]/n;
if (Amax<a[j]) Amax=a[j];
if (Bmax<b[j]) Bmax=b[j];
if (Amin>a[j]) Amin=a[j];
if (Bmin>b[j]) Bmin=b[j];
}
Amax=fabs(Amax);
Amin=fabs(Amin);
Bmax=fabs(Bmax);
Bmin=fabs(Bmin);
if(Amax>Amin) AM=Amax;
else AM=Amin;
if(Bmax>Bmin) BM=Bmax;
else BM=Bmin;
if(AM>BM) Pk=AM;
else Pk=BM;
Pp=-Pk;
double q=0;
q=640/n;
outtextxy(40,20,"a[k]");
line (10,240,630,240);
line (10,10,10,470);
for (j=0;j<n;j++)
{
ag[j]=Pk-(a[j]-Pp)*(Pgk-Pgp)/(Pk-Pp);
G[j]=15+j*q;
line (G[j],240,G[j],ag[j]);
}
getch();
cleardevice();
outtextxy(40,20,"b[k]");
line (10,240,630,240);
line (10,10,10,470);
for (j=0;j<n;j++)
{
bg[j]=Pk-(b[j]-Pp)*(Pgk-Pgp)/(Pk-Pp);
G[j]=15+j*q;
line (G[j],240,G[j],bg[j]);
}
getch();
cleardevice();
double Gser=0, Sigma=0;
for(j=0;j<n;j++)
{
Gser=Gser+sqrt(a[j]*a[j]+b[j]*b[j]);
}
Gser=Gser/n;
Sigma=Gser/1.25;
double pGa=0,dG,Ga=0;
dG=Ga/n;
i=0;
for(Ga=0;Ga<0.1;Ga+=0.1/n)
{
Y[i]=(Ga*exp(-(pow(Ga,2))/(2*pow(Sigma,2))))/(pow(Sigma,2));
X[i]=i;
i++;
}
dG=1000/n;
while(pGa>=Ga)
{
pGa=pGa+(float)Ga*exp(-(Ga*Ga)/(2*Sigma*Sigma))/(Sigma*Sigma)*dG;
Ga+=dG;
}
float beta=Ga/Gser;
for(j=0;j<n;j++)
{
if(sqrt(a[j]*a[j]+b[j]*b[j])>=beta*Ga)
{
a[j]=Gser/sqrt(2);
b[j]=Gser/sqrt(2);
}
X[j]=j;
}
Xp=0; Yp=0; Xk=0; Yk=0;
Xgp=20, Xgk=620, Ygp=20, Ygk=460;
line(Xgp,Ygk,Xgk,Ygk);
line(Xgp,Ygp,Xgp,Ygk);
outtextxy(40,20,"DataNois after filtering");
for(i=0;i<n;i++)
{
Y[i]=a[0]/2;
for(j=1;j<n;j++)
{
Y[i]=Y[i]+a[j]*cos(2*M_PI*j*i/n)+sin(2*M_PI*j*i/n);
X[i]=0+i*h;
}
if (X[i]>Xk) Xk=X[i];
if (X[i]<Xp) Xp=X[i];
if (Y[i]>Yk) Yk=Y[i];
if (Y[i]<Yp) Yp=Y[i];
}
for (i=0;i<n;i++)
{
Xg=Xgp+(X[i]-Xp)*(Xgk-Xgp)/(Xk-Xp);
Yg=Ygk-(Y[i]-Yp)*(Ygk-Ygp)/(Yk-Yp);
putpixel (Xg,Yg,2);
}
getch();
cleardevice();
closegraph();
return 0;
}
Результати роботи програми:
Рис.1 Графічне представлення вимірювального сигналу.
Рис.2 Графічне представлення завади.
Рис.3 Графічне представлення вимірювального сигналу та завади.
Рис.4 Графічне представлення спектру сигналу.
Рис.5 Графічне представлення вимірювального сигналу після фільтрації.