# A Fast Power Grid Frequency Estimation Approach Using Frequency-Shift Filtering using numerical example

I have read the paper A Fast Power Grid Frequency Estimation Approach Using Frequency Shift Filtering.

I want to prove it using a numerical example, I have problem how to calculate the hp(n) values,

Let's suppose that I have a sample signal which I have sampled with an ADC and for simplicity I will use Matlab to make the signal. (the signal has a 50.1HZ frequency) and I want to estimate it with this paper approach,

fs = 12800;                    % Sampling frequency (samples per second)
dt = 1/fs;                   % seconds per sample
StopTime = 0.02;             % 20 mili seconds or a 1 full cycle of 50Hz AC
t = (0:dt:StopTime-dt)';     % seconds
F = 50.1;                    % Sine wave frequency (hertz)
X = sin(2*pi*F*t);           % Generate sine wave


Now How should I calculate the hp(n) values using MATLAB? Finally I should convolve the result with this sample sin wave to get 50.1 Hz estimation.

• The article seems to be behind a paywall... – Ben Feb 16 at 12:45
• From the look of it, they seem to perform the Park transform with a fixed-frequency (50 Hz) in your case. Since, the frequency is not exactly 50 Hz, the d-q signals won't be DC but will have a frequency of f - 50 Hz or 0.1 Hz in your case. – Ben Feb 16 at 12:51
• Thanks, then How should I generate the H and do convolution? – ASiDesigner Feb 16 at 14:33
• you can check out the Paper in here sci-hub.do/10.1109/tpwrs.2019.2892599 – ASiDesigner Feb 16 at 15:00
• The filter is a moving-average, h[n] is explained on page 2 $h_{av}$ – Ben Feb 16 at 16:54

Basically they cascade moving-average filters, or CAF.

a Moving-average filter has the following coefficients

$$h_{av}(m) = 1/M, m = 0, 1.. M-1$$

In Matlab it would look like this

h = ones(1,M)./M;
y = filter(h,1,x);


If you want to cascade three moving-average filters

h = ones(1,M)./M;
y1 = filter(h,1,x);
y2 = filter(h,1,y1);
y3 = filter(h,1,y2);


M should be equal to the number of samples per period. So if your frequency is approximately 50 Hz and your sampling frequency is 12800 Hz then M = 256;.

• Thanks, for the info, so how hp is calculated, so we do the convolution and estimate that this example X signal has 50.1Hz frequency? – ASiDesigner Feb 18 at 11:16
• No, the rest of the processing required is explained in the paper. You need to shift the frequency before filtering. – Ben Feb 18 at 12:34
• Dear Ben thanks for the feedback, would you please tell me How to do it in matlab? I'm a bit confused how to do it. thanks – ASiDesigner Feb 20 at 6:54
• How should I do it in matlab? – ASiDesigner Feb 23 at 7:55
• Which part specifically? – Ben Feb 23 at 14:21

So, Now I have done a little bit more, Now I'm generating a complex reference sine wave, The things that I have done are as follow

clear
clc
%ADC sin wave, it represent the data that I have captured with the ADC
%and want to estimate the frequency
grid on
hold on

%ref sin wave with 50Hz complex value
Refsin = dsp.SineWave(1,50);
Refsin.ComplexOutput = true;
Refsin.SampleRate  = 12800;
Refsin.SamplesPerFrame = 256;
ref = Refsin();
plot(imag(ref))

%calculating  xs
xs = ref .* ADC;

%calculate the CAF
M = 256;
h = ones(1,M)./M;
y1 = filter(h,1,ref);
y2 = filter(h,1,y1);
y3 = filter(h,1,y2);

%doing the actual final convolution
F = conv2(xs, y3);

%calculating the final result
Fest = 50 - (angle(F(4)) - angle(F(3)) * 256*50)/(2*pi*90)


the answer is clearly wrong! so what I have done wrong?