procedure NOTCH(var WorkArr1: diarrWork;
iDuration, iRate: integer; realF0, realdF: real; iOrd: integer);
{
===============================================================================
PURPOSE : Notch frequency filtering of data with Butterword
filter of the second order.
I/O data : Input: WorkArr1 - array with initial data;
iDuration - length of array;
iRate - rate of DAQ;
realF0 - cut frequency;
realDf - frequency band;
iOrd - number of averaging iterations.
Output: WorkArr1 -array with filtered data.
-------------------------------------------------------------------------------
BEFORE : Array WorkArr1 with initial data.
AFTER : Array Workarr1 with filtered data.
}
type
arrealCoef = array [1..50] of real;
darrealPoint = ^arrealPoint;
arrealPoint = array [1..cNPointMax] of real;
var
arA, arB: darrealPoint;
realDT, realBZero, realB1, realB2: real;
arA1, arA2: arrealCoef;
NRepeat, i, i3, NPoint: integer;
(*
type arcoef = array [1..50] of single;
arfreq = array [1..50] of single;
arsingle = array [0..DeltMax-1] of single;
ard = ^arsingle;
var
aaa,bbb:ard;
dt:single;
a1,a2:arcoef;
bzero,b1,b2:single;
nff:integer;
i,i3:integer;
*)
function atan2(ar, ai: real): real;
{
===============================================================================
PURPOSE : Calculation of ArcTan with the main value corresponding
to signs of numerator and denominator.
I/O data : Input: ar - numerator;
ai - denominator.
Output: function atan2 with main value.
-------------------------------------------------------------------------------
BEFORE :
AFTER :
}
begin
if ar = 0 then
begin
if ai > 0 then atan2 := 0 else atan2 := -Pi;
end
else begin { ar <> 0 }
if ai = 0 then
begin
if ar > 0 then atan2 := pi / 2.0
else atan2 := -pi / 2.0;
end { ai = 0 }
else begin { ai <> 0 }
atan2 := arctan(ar / ai);
end; { ai <> 0 }
end; { ar <> 0 }
end; { function atan2 }
procedure ttran(var arA1, arA2, arB0, arB1, arB2: arrealCoef; m, npo: integer;
realT: real; var freq, abz, phs: arrealCoef);
{
===============================================================================
PURPOSE : Solution of charateristic equation for
given f0 and df.
I/O data : Input: arA1,arA2,arB0,arB1,arB2 - arrays with filter coefficients
realT - normalized time step between DAQ readings;
m - order of filter;
npo = 1 - fixed value of iterations;
freq array with arguments;
Output: abz - array with solutions' amplitudes;
phs - array with solutions' phases.
-------------------------------------------------------------------------------
BEFORE :
AFTER : Calculated solutions of charateristic equation.
}
var
fact, add, prev, fd, s1, c1, a, s2, c2, absa, phsa,
ar, ai, anm, pnd, dum, ang, cur, test: real;
ip, i, j: integer;
begin
fact := 2 * Pi * realT;
ip := m - m div 2;
add := 0;
prev := 0;
for i := 1 to npo do
begin
fd := freq[i] * fact;
s1 := sin(fd);
c1 := cos(fd);
a := 2.0 * fd;
s2 := sin(a);
c2 := cos(a);
absa := 1.0;
phsa := 0;
for j := 1 to ip do
begin
ar := arB0[j] + arB1[j] * c1 + arB2[j] * c2;
ai := -arB1[j] * s1 - arB2[j] * s2;
anm := ar * ar + ai * ai;
pnd := 0;
if (ai <> 0) or (ar <> 0) then pnd := atan2(ai, ar);
ar := 1.0 + arA1[j] * c1 + arA2[j] * c2;
ai := -arA1[j] * s1 - arA2[j] * s2;
absa := absa * anm / (ar * ar + ai * ai);
dum := 0;
if (ai <> 0) or (ar <> 0) then dum := atan2(ai, ar);
phsa := phsa + pnd - dum;
end; { for j := 1 to ip do }
abz[i] := 10.0 * ln(absa + 1.0e-30) / ln(10.0);
ang := phsa * 180.0 / Pi;
cur := ang + add;
test := cur - prev;
if abs(test) >= 180 then
begin
if (test < 0) then
begin
add := add + 360;
cur := cur - 360;
end { if (test < 0) then }
else begin
add := add - 360;
cur := cur + 360;
end; { if (test >= 0) then }
end; { if abs(test) >= 180 then }
prev := cur;
phs[i] := cur;
end; { for i := 1 to npo do }
Exit;
end; { procedure ttran(var arA1, arA2, arB0, arB1, arB2: arrealCoef; m, npo: integer;
realT: real; var freq, abz, phs: arrealCoef); }
procedure Bnps(mm: integer; var t, bw, fc: real;
var arA1, arA2: arrealCoef; var bzero, b10, b20: real);
{
===============================================================================
PURPOSE : Calculation of Notch filter parameters for
given f0 and df.
I/O data : Input: mm - order of filter;
t - time step between DAQ readings;
bw - cut frequency;
fc - frequency band.
Output: arA1, arA2 - arrays with filter coefficients;
bzero,b10,b20 - filter coefficient.
-------------------------------------------------------------------------------
BEFORE :
AFTER : Calculated filter coefficients.
}
var
arB0, arB1, arB2, freq, abz, phs: arrealCoef;
a, b, c, d, e, g, h, fact, wedge, sector, cc, ss, ang2, fn,
ang, htran: real;
i, m, m1: integer;
begin
arB0[1] := 1.0;
arB1[1] := 0;
arB2[1] := 1.0;
fact := 2.0 * pi * t * bw;
ang2 := 2.0 * pi * t * fc;
cc := cos(ang2) * cos(fact);
arB1[1] := -2.0 * cc;
freq[1] := 0;
ss := sin(ang2) * sin(fact);
m := mm;
m1 := m div 2;
a := m;
htran := 0;
sector := 2 * pi / a;
wedge := sector / 2.0;
for i := 1 to m1 do
begin
fn := i - 1;
ang := fn * sector + wedge;
a := ss * cos(ang) + cc;
b := ss * sin(ang);
c := 1.0 - (a * a + b * b);
d := 0.5 * (-c + sqrt(c * c + 4.0 * b * b));
e := sqrt(d + 1.0) + sqrt(d);
if (d = 0) then //VV01January2014
g := 2.0 * sqrt(1.0 - c) / e //VV01January2014
else //VV01January2014
g := 2.0 * sqrt(1.0 - b * b / d) / e;
if a < 0 then g := -g;
h := -1.0 / (e * e);
arA1[i] := -g;
arA2[i] := -h;
ttran(arA1, arA2, arB0, arB1, arB2, 2, 1, 1.0, freq, abz, phs);
htran := htran + abz[1];
end; { for i := 1 to m1 do }
bzero := exp(ln(10.0) * (-htran / (20.0 * m1)));
b10 := arB1[1] * bzero;
b20 := arB2[1] * bzero;
exit;
end; { procedure Bnps }
begin {procedure NOTCH(var WorkArr1: diarrWork;
iDuration, iRate: integer; realF0, realdF: real; iOrd: integer); }
NPoint := iRate * iDuration;
new(arA);
new(arB);
realDt := 1.0 / iRate;
// realDt := 0.001;
Bnps(2, realDt, realDf, realF0, arA1, arA2, realBzero, realB1, realB2);
// Lpsb(2, realDt, realF0, arA1, arA2, realBzero); { calculation
// of filter coefficients }
for i := 1 to NPoint do arA^[i] := WorkArr1^[i];
NRepeat := (iOrd div 2) * 2;
if NRepeat <= 0 then NRepeat := 2;
for i3 := 1 to NRepeat do
begin { cycle by iterations }
//VV01January2014 arB^[1] := (realBzero * arA^[1]);
//VV01January2014 arB^[2] := realBzero * arA^[2] - arA1[1] * arB^[1] + realB1 * arA^[1];
arB^[1] := (arA^[1]); //VV01January2014
arB^[2] := (arA^[2]); //VV01January2014
{--- filtering ---}
for i := 3 to NPoint do
arB^[i] := realBzero * (arA^[i] + arA^[i - 2]) - arA1[1] * arB^[i - 1]
- arA2[1] * arB^[i - 2] + realB1 * arA^[i - 1];
for i := 1 to NPoint do arA^[NPoint - i + 1] := arB^[i]; {inverse sequence}
end; { for i3 := 1 to NRepeat do }
i3 := Round(iRate / Realf0);
for i := i3 downto 1 do
arB^[i] := arB^[i + 1];
for i := i3 to NPoint do
arB^[i] := arB^[i - 1];
(*
new(aaa);
new(bbb);
dt:=1.0/fdiscr;
Bnps(2,dt,df,f0,a1,a2,bzero,b1,b2); // filter coefficients calculation
for i:=0 to deltk-1 do aaa^[i]:=WorkArr1^[i];
if order = 1 then nff:= 1 else
nff:=(order div 2)*2;
for i3:=1 to nff do begin // cycle by iterations
bbb^[0]:=(bzero*aaa^[0]);
bbb^[1]:=bzero*aaa^[1]-a1[1]*bbb^[0]+b1*aaa^[0];
for i:=2 to deltk-1 do
bbb^[i]:=bzero*(aaa^[i]+aaa^[i-2])-a1[1]*bbb^[i-1]-
a2[1]*bbb^[i-2]+b1*aaa^[i-1];
if nff = 1 then
for i:=0 to deltk-1 do aaa^[i]:=bbb^[i]
else
for i:=0 to deltk-1 do aaa^[deltk-i-1]:=bbb^[i]; // inverse sequency
end; {--i3--}
for i:=0 to deltk-1 do begin
if (aaa^[i] >= 32766.0) then aaa^[i]:=32766.0
else if aaa^[i] < -32767.0 then aaa^[i]:=-32767.0;
WorkArr1^[i]:=round(aaa^[i]);
end; // i
dispose(aaa);
dispose(bbb);
*)
for i := 1 to NPoint do
WorkArr1^[i] := Round(arA^[i]);
Dispose(arA);
Dispose(arB);
end;{ procedure NOTCH }
最新发布
文章探讨了非侵入性负载监测（NILM）问题，将其视为时间和序列分析的优化问题。通过使用seq2seq模型和LSTM网络，然后提出seq2point模型来改进预测准确性。seq2point模型利用1D卷积神经网络集中处理滑动窗口的中点，从而提供更精确的设备功率预测。损失函数用于优化模型，而模型结构包括一系列1D卷积层和全连接层。
最低0.47元/天 解锁文章
扫一扫
8
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
