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

This commit is contained in:
2026-06-21 00:52:59 +03:00
parent 81a3e17f7e
commit 6f89664522
10 changed files with 286 additions and 99041 deletions
File diff suppressed because it is too large Load Diff
@@ -1,398 +0,0 @@
// Статические переменные
funcprot(0);
Fx = 330000; // Измеряемая частота
F1 = Fx*1.1; // Максимальная частота
F11 = Fx*0.9; // Минимальная чстота
A1 = 0.25; // Амплитуда измеряемой величины
Fp = 100; // Частота помехи
Ap = 1; // Амплитуда помехи
// Tm = 10e-3; // Время выборки
Fd = 168/9*1e6; // Частота дискретизации
NADC = 8; // Количество бит АЦП
nmem = 2^16; // количество ячеек памяти
// определим время выборки
Td = 1/Fd;
Tm = nmem * Td
nn=int(Tm/Td);
Aop=(1.25+1);
Amax=(1.25+1)+Aop;
AmaxADC = Amax*1.2;
StepADC = AmaxADC / (2^NADC);
// Fi1 - начальная фаза исходного сигнала
// Fip - начальная фаза помехи
// Генерация выборок исходного сигнала
Fi1=0
function s1 = ADC1()
// Fi1 = rand()*2*%pi;
//Fi1=Fi1+%pi/100;
disp(Fi1);
Fip = rand()*2*%pi;
Omega1 = 2*%pi*F1;
Omegap = 2*%pi*Fp;
s1=zeros(nn);
// Отклонение частоты генератора
if rand()>0.5 then
Td123 = Td*1.0001;
else
Td123 = Td*0.9999;
end
disp(1/Td123/1e6)
for ii=1:nn
ti = Td123*ii
sTMP = A1*cos(Omega1*ti+Fi1) + Ap*cos(Omegap*ti+Fip) + Aop;
s1(ii) = int(sTMP/StepADC)*StepADC; // Округление АЦП
end
endfunction
// Формирование целочисленного массива 256 уровней
function s1 = ADC2()
// Fi1 = rand()*2*%pi;
//Fi1=Fi1+%pi/100;
disp(Fi1);
Fip = rand()*2*%pi;
Omega1 = 2*%pi*F1;
Omegap = 2*%pi*Fp;
s1=zeros(nn);
// Отклонение частоты генератора
if rand()>0.5 then
Td123 = Td*1.0001;
else
Td123 = Td*0.9999;
end
disp(1/Td123/1e6)
for ii=1:nn
ti = Td123*ii
sTMP = A1*cos(Omega1*ti+Fi1) + Ap*cos(Omegap*ti+Fip) + Aop;
s1(ii) = int(sTMP/StepADC); // Округление АЦП
end
endfunction
// Медианный фильтр по пяти точкам
function s2 = Median5(s1)
nn=size(s1,1);
for ii = 1:(nn-5)
v1(1)=s1(ii);
v1(2)=s1(ii+1);
v1(3)=s1(ii+2);
v1(4)=s1(ii+3);
v1(5)=s1(ii+4);
v2 = gsort(v1);
s2(ii) = v2(3);
end
s2(nn-4)=s1(nn-4);
s2(nn-3)=s1(nn-3);
s2(nn-2)=s1(nn-2);
s2(nn-1)=s1(nn-1);
s2(nn)=s1(nn);
endfunction
// Медианный фильтр по трём точкам
function s2 = Median3(s1)
nn = size(s1,1);
for ii = 1:(nn-5)
v1(1)=s1(ii);
v1(2)=s1(ii+1);
v1(3)=s1(ii+2);
v1(4)=s1(ii+3);
v1(5)=s1(ii+4);
v2 = gsort(v1);
s2(ii) = v2(3);
end
s2(nn-4)=s1(nn-4);
s2(nn-3)=s1(nn-3);
s2(nn-2)=s1(nn-2);
s2(nn-1)=s1(nn-1);
s2(nn)=s1(nn);
endfunction
// Скользящее среднее по пяти точкам
function s2 = Mov_avg(s1)
nn = size(s1,1);
s2 = zeros(1:nn);
x1 = 0;
for i=1:5
x1 = x1 + s1(i);
s2(i) = s1(i);
end
for i=5:nn-1
s2(i) = int(x1/5);
x1 = x1-s1(i-4);
x1 = x1+s1(i+1);
end
s2(nn) = s1(nn);
endfunction
// Средний размах + Оценка периода
function [Fest, Sum1, s3] = Filter2(s2)
nx = int((2/F11)/Td);
nx2 = int((1/4/F11)/Td);
// disp(nx);
// disp(nx2);
if nx2<2 then
nx2 = 2;
end
s5 = zeros(nn);
// for ii=1:nn-nx+1
// s4(1:nx) = s2(ii:ii+nx-1)
// vmax=max(s4);
// vmin=min(s4);
// vref=(vmax+vmin)/2;
// s5(ii) = s2(ii) - vref;
// s3(ii) = vmax-vmin;
// end
for ii=0:int(nn/nx)-1
s4(1:nx) = s2(1+ii*nx:ii*nx+nx)
vmax=max(s4);
vmin=min(s4);
vref=(vmax+vmin)/2;
s5(1+ii*nx:ii*nx+nx) = s2(1+ii*nx:ii*nx+nx) - vref;
s3(1+ii*nx:ii*nx+nx) = vmax-vmin;
end
Sum1=0;
for ii=1:nn-nx;
Sum1=Sum1+s3(ii);
end
Sum1=Sum1/(nn-nx);
// Считаем периоды и оцениваем частоту
nper = 0;
jj=1;
while s5(jj)>0
jj=jj+1
end
while s5(jj)<0
jj=jj+1
end
kk1 = jj
while jj<nn-nx
jj = jj + nx2;
while (jj<nn-nx) && (s5(jj)<0)
jj=jj+1;
end
nper = nper +1;
kk2 = jj;
jj = jj + nx2;
while (jj<nn-nx) &&(s5(jj)>0)
jj=jj+1;
end
end
Fest = 1/(((kk2-kk1)*Td)/nper);
endfunction
function gg=MyHist1(sw1)
// Формируем гистограмму из массива
// знаем, что разрядность АЦП 8 бит - 256 значений
gg=zeros(1:256);
for i=1:size(sw1,1)
nn=int(sw1(i)/StepADC);
gg(nn)=gg(nn)+1;
end
endfunction
function gg=MyHist2(sw1)
// Формируем гистограмму из массива
// знаем, что разрядность АЦП 8 бит - 256 значений
gg=zeros(1:256)
for i=1:size(sw1,1)
gg(sw1(i))=gg(sw1(i))+1
end
endfunction
function [gg, Fest, Vest]=Estim2(s1)
// !!!!!!! Обработка целочисленных данных !!!!!!!!
// Размер окна примерно 10 периодов на минимальной частоте
wn = int(Fd/F11)*10;
// Объём входных данных
nmax = size(s1,1);
ymin = 256;
ymax = 0;
ypp = 0;
yhist = zeros(1:256);
xhist = zeros(1:256);
prlag = 8; // Защитный интервал
prtmp = 8;
Tr2 = 0;
xzero = 0;
xn1 = 1;
xn2 = wn;
// Первый цикл - инициализация
for i=xn1:xn2
if ymin>s1(i) then ymin = s1(i) end;
if ymax<s1(i) then ymax = s1(i) end;
end
ypp = ymax-ymin;
yhist(ypp) = yhist(ypp)+1;
trh = int(ypp/2)+ymin;
// Первый перход через 0
Tr3 = 0;
for i=xn1:xn2
// Поиск переходов через 0
stmp = s1(i);
if prtmp >= prlag then
if Tr2 == 1 then
if stmp>trh then
prtmp = 0;
if Tr3 == 1 then
xhist(i-xzero)=xhist(i-xzero)+1;
end
xzero = i;
Tr3 = 1;
Tr2 = 0;
end
else
if stmp<trh then
prtmp = 0;
Tr2 = 1;
end
end
end
prtmp = prtmp+1;
end;
while xn2<nmax
xn1 = xn1+1;
xn2 = xn2+1;
stmp = s1(xn2);
Tr1 = 0;
if s1(xn2)>ymax then
ymax = stmp;
elseif s1(xn1-1)==ymax then
Tr1=1;
end
if s1(xn2)<ymin then
ymin = stmp;
elseif s1(xn1-1)==ymin then
Tr1 = 1;
end
if Tr1==1 then
ymin = 256;
ymax = 0;
for i=xn1:xn2
if ymin>s1(i) then ymin = s1(i) end;
if ymax<s1(i) then ymax = s1(i) end;
end
end
ypp = ymax-ymin;
yhist(ypp) = yhist(ypp)+1;
// Поиск переходов через 0
trh = int(ypp/2)+ymin;
if prtmp >= prlag then
if Tr2 == 1 then
if stmp>trh then
prtmp = 0;
xhist(xn2-xzero)=xhist(xn2-xzero)+1;
xzero = xn2;
Tr2 = 0;
end
else
if stmp<trh then
prtmp = 0;
Tr2 = 1;
end
end
end
prtmp = prtmp+1;
end
gg = yhist;
// Оценка гистограмм
x=0;
y=0;
for i=1:256
if yhist(i)>y then
y=yhist(i);
Vest=i;
end
if xhist(i)>x then
x=xhist(i);
Fest=i;
end
end
endfunction
rand("uniform");
rand("seed");
//Fi1 = 0.1570796-%pi/1000;
//Fi1 = 0
//for rkk=1:100
// disp(rkk);
// Fi1=Fi1+%pi/1000;
// s1 = ADC1();
// s2 = Filter1(s1);
// [Fest, Sum1, s3] = Filter2(s2);
// disp((Sum1/2-A1)/A1*100);
// disp((Fest-F1)/F1*100);
//end
//s1 = ADC1();
//sw1 = s1(1:1024);
//gg = MyHist1(sw1);
//plot2d(gg);
fl1(1)="c:\TestADC\10khz_809mv_without_R.txt";
fl1(2)="c:\TestADC\10khz_809mv_with_R.txt";
fl1(3)="c:\TestADC\300khz_201mv_without_R.txt";
fl1(4)="c:\TestADC\300khz_201mv_with_R.txt";
fl1(5)="c:\TestADC\300khz_802mv_without_R.txt";
fl1(6)="c:\TestADC\300khz_802mv_with_R.txt";
fl1(7)="c:\TestADC\400khz_201mv_without_R.txt";
fl1(8)="c:\TestADC\400khz_201mv_with_R.txt";
fl1(9)="c:\TestADC\400khz_802mv_without_R.txt";
fl1(10)="c:\TestADC\400khz_802mv_with_R.txt";
fl1(11)="c:\TestADC\500khz_201mv_without_R.txt";
fl1(12)="c:\TestADC\500khz_201mv_with_R.txt";
fl1(13)="c:\TestADC\500khz_802mv_without_R.txt";
fl1(14)="c:\TestADC\500khz_802mv_with_R.txt";
uinp(1) = 809;
uinp(2) = 809;
uinp(3) = 201;
uinp(4) = 201;
uinp(5) = 802;
uinp(6) = 802;
uinp(7) = 201;
uinp(8) = 201;
uinp(9) = 802;
uinp(10) = 802;
uinp(11) = 201;
uinp(12) = 201;
uinp(13) = 802;
uinp(14) = 802;
for fff=3:14
s1 = read(fl1(fff),-1,1);
s2 = Mov_avg(s1);
h=scf(1);
clf;
plot2d(s1);
[gg, Fest, Vest]=Estim2(s2);
Err1 = Vest/255*1000;
Err1 = (Err1-uinp(fff))/uinp(fff)*100;
disp(fl1(fff), Err1);
h=scf(2);
clf;
plot2d(gg);
pause;
end
//sw1 = s1(1:1500);
//gg = MyHist2(sw1);
//plot2d(s1);
//s1 = ADC2();
@@ -1,339 +0,0 @@
// Статические переменные
funcprot(0);
Fd = 168/9*1e6; // Частота дискретизации
NADC = 8; // Количество бит АЦП
nmem1 = 2^16; // количество ячеек памяти
nmem2 = nmem1 -5; // после медианного фильтра
// определим время выборки
Td = 1/Fd;
Tm = nmem * Td
// Медианный фильтр по пяти точкам
function s2 = Median5(s1, nn)
s2 = zeros(1:nn-5); // Обнуление выходного массива
for ii = 1:(nn-5)
v1(1)=s1(ii);
v1(2)=s1(ii+1);
v1(3)=s1(ii+2);
v1(4)=s1(ii+3);
v1(5)=s1(ii+4);
v2 = gsort(v1);
s2(ii) = v2(3);
end
endfunction
// Скользящее среднее по пяти точкам
function s2 = Mov_avg(s1, nn)
s2 = zeros(1:nn); // Обнуление выходного массива
x1 = 0;
for ii = 1:5
x1 = x1 + s1(ii);
s2(ii) = s1(ii);
end
for ii = 1:(nn-5)
s2(ii) = int(x1/5);
x1 = x1 - s1(ii+1);
x1 = x1 + s1(ii+2);
end
endfunction
// Получение оценки максимума и минимума по гистограмме
function [ymin, ymax] = Histo1(s2)
// Формируем гистограмму из массива
// знаем, что разрядность АЦП 8 бит - 256 значений
// тип данных gg unsigned short
gg = zeros(1:256); // Обнуление массива
for ii = 1:nmem2
gg(s2(ii)) = gg(s2(ii)) + 1;
end
// Поиск максимумов в гистограмме
ymin = 0;
ymax = 0;
y = 0;
for ii = 1:124 // Серединку не обрабатываем
if gg(ii) > y then
y = gg(ii);
ymin = ii;
end
y = 0;
for ii = 130:256
if gg(ii) > y then
y = g(ii);
ymax = ii;
end
endfunction
// Последовательный поиск максимумов и минимумов, оценка пиков
// и вычисление среднего размаха
// Средний размах + Оценка периода
function [Fest, Sum1, s3] = Filter2(s2)
nx = int((2/F11)/Td);
nx2 = int((1/4/F11)/Td);
// disp(nx);
// disp(nx2);
if nx2<2 then
nx2 = 2;
end
s5 = zeros(nn);
// for ii=1:nn-nx+1
// s4(1:nx) = s2(ii:ii+nx-1)
// vmax=max(s4);
// vmin=min(s4);
// vref=(vmax+vmin)/2;
// s5(ii) = s2(ii) - vref;
// s3(ii) = vmax-vmin;
// end
for ii=0:int(nn/nx)-1
s4(1:nx) = s2(1+ii*nx:ii*nx+nx)
vmax=max(s4);
vmin=min(s4);
vref=(vmax+vmin)/2;
s5(1+ii*nx:ii*nx+nx) = s2(1+ii*nx:ii*nx+nx) - vref;
s3(1+ii*nx:ii*nx+nx) = vmax-vmin;
end
Sum1=0;
for ii=1:nn-nx;
Sum1=Sum1+s3(ii);
end
Sum1=Sum1/(nn-nx);
// Считаем периоды и оцениваем частоту
nper = 0;
jj=1;
while s5(jj)>0
jj=jj+1
end
while s5(jj)<0
jj=jj+1
end
kk1 = jj
while jj<nn-nx
jj = jj + nx2;
while (jj<nn-nx) && (s5(jj)<0)
jj=jj+1;
end
nper = nper +1;
kk2 = jj;
jj = jj + nx2;
while (jj<nn-nx) &&(s5(jj)>0)
jj=jj+1;
end
end
Fest = 1/(((kk2-kk1)*Td)/nper);
endfunction
function [gg, Fest, Vest]=Estim2(s1)
// !!!!!!! Обработка целочисленных данных !!!!!!!!
// Размер окна примерно 10 периодов на минимальной частоте
wn = int(Fd/F11)*10;
// Объём входных данных
nmax = size(s1,1);
ymin = 256;
ymax = 0;
ypp = 0;
yhist = zeros(1:256);
xhist = zeros(1:256);
prlag = 8; // Защитный интервал
prtmp = 8;
Tr2 = 0;
xzero = 0;
xn1 = 1;
xn2 = wn;
// Первый цикл - инициализация
for i=xn1:xn2
if ymin>s1(i) then ymin = s1(i) end;
if ymax<s1(i) then ymax = s1(i) end;
end
ypp = ymax-ymin;
yhist(ypp) = yhist(ypp)+1;
trh = int(ypp/2)+ymin;
// Первый перход через 0
Tr3 = 0;
for i=xn1:xn2
// Поиск переходов через 0
stmp = s1(i);
if prtmp >= prlag then
if Tr2 == 1 then
if stmp>trh then
prtmp = 0;
if Tr3 == 1 then
xhist(i-xzero)=xhist(i-xzero)+1;
end
xzero = i;
Tr3 = 1;
Tr2 = 0;
end
else
if stmp<trh then
prtmp = 0;
Tr2 = 1;
end
end
end
prtmp = prtmp+1;
end;
while xn2<nmax
xn1 = xn1+1;
xn2 = xn2+1;
stmp = s1(xn2);
Tr1 = 0;
if s1(xn2)>ymax then
ymax = stmp;
elseif s1(xn1-1)==ymax then
Tr1=1;
end
if s1(xn2)<ymin then
ymin = stmp;
elseif s1(xn1-1)==ymin then
Tr1 = 1;
end
if Tr1==1 then
ymin = 256;
ymax = 0;
for i=xn1:xn2
if ymin>s1(i) then ymin = s1(i) end;
if ymax<s1(i) then ymax = s1(i) end;
end
end
ypp = ymax-ymin;
yhist(ypp) = yhist(ypp)+1;
// Поиск переходов через 0
trh = int(ypp/2)+ymin;
if prtmp >= prlag then
if Tr2 == 1 then
if stmp>trh then
prtmp = 0;
xhist(xn2-xzero)=xhist(xn2-xzero)+1;
xzero = xn2;
Tr2 = 0;
end
else
if stmp<trh then
prtmp = 0;
Tr2 = 1;
end
end
end
prtmp = prtmp+1;
end
gg = yhist;
// Оценка гистограмм
x=0;
y=0;
for i=1:256
if yhist(i)>y then
y=yhist(i);
Vest=i;
end
if xhist(i)>x then
x=xhist(i);
Fest=i;
end
end
endfunction
rand("uniform");
rand("seed");
//Fi1 = 0.1570796-%pi/1000;
//Fi1 = 0
//for rkk=1:100
// disp(rkk);
// Fi1=Fi1+%pi/1000;
// s1 = ADC1();
// s2 = Filter1(s1);
// [Fest, Sum1, s3] = Filter2(s2);
// disp((Sum1/2-A1)/A1*100);
// disp((Fest-F1)/F1*100);
//end
//s1 = ADC1();
//sw1 = s1(1:1024);
//gg = MyHist1(sw1);
//plot2d(gg);
fl1(1)="c:\TestADC\10khz_809mv_without_R.txt";
fl1(2)="c:\TestADC\10khz_809mv_with_R.txt";
fl1(3)="c:\TestADC\300khz_201mv_without_R.txt";
fl1(4)="c:\TestADC\300khz_201mv_with_R.txt";
fl1(5)="c:\TestADC\300khz_802mv_without_R.txt";
fl1(6)="c:\TestADC\300khz_802mv_with_R.txt";
fl1(7)="c:\TestADC\400khz_201mv_without_R.txt";
fl1(8)="c:\TestADC\400khz_201mv_with_R.txt";
fl1(9)="c:\TestADC\400khz_802mv_without_R.txt";
fl1(10)="c:\TestADC\400khz_802mv_with_R.txt";
fl1(11)="c:\TestADC\500khz_201mv_without_R.txt";
fl1(12)="c:\TestADC\500khz_201mv_with_R.txt";
fl1(13)="c:\TestADC\500khz_802mv_without_R.txt";
fl1(14)="c:\TestADC\500khz_802mv_with_R.txt";
uinp(1) = 809;
uinp(2) = 809;
uinp(3) = 201;
uinp(4) = 201;
uinp(5) = 802;
uinp(6) = 802;
uinp(7) = 201;
uinp(8) = 201;
uinp(9) = 802;
uinp(10) = 802;
uinp(11) = 201;
uinp(12) = 201;
uinp(13) = 802;
uinp(14) = 802;
for fff=3:14
s1 = read(fl1(fff),-1,1);
s2 = Mov_avg(s1);
h=scf(1);
clf;
plot2d(s2);
[gg, Fest, Vest]=Estim2(s2);
Err1 = Vest/255*1000;
Err1 = (Err1-uinp(fff))/uinp(fff)*100;
disp(fl1(fff), Err1);
h=scf(2);
clf;
plot2d(gg);
pause;
end
//sw1 = s1(1:1500);
//gg = MyHist2(sw1);
//plot2d(s1);
//s1 = ADC2();