Recovering shape of periodically sampled signal in near resonance conditions

Periodic signal with frequency Fs has sampled with frequency Fo. How can we recover the signal's shape if Fs is close to Fo or one of its harmonics, so the period of a beating frequency takes hundreds of thousands of samples?

Thanks for detailed answer

[Sulaymon Eshkabilov]

This is a good theoretical explantion.

[Valeriy]

Thank you very much, @Sulaymon Eshkabilov, this is really good article.

Do you have information about Matlab realization of such algorithms?

[Sulaymon Eshkabilov]

No, I don't. By following their guidelines, one can write a decent code to get the job done.

Answer 1

[Valeriy]

@Sulaymon Eshkabilov, I'm not working in mathematics area, so I'm not familiar with common abbreviations and I have a difficulty to apply ideas of cited publication to my case. Please, be so kind to help me.

From the cited publication:

At first, in my case there are no frequency modulation, even more, signal frequency is constant and I know its value with good precision. So, if I understand correctly, theta dependence in my case is absent.

Second, sampling of my signal realized with constant, uniform spatial frequency, so how to understand interpolation of the original signal to uniform grid? It is already in uniform grid.

Is shape function from this publication shape of my periodic signal?

And, what does it mean, interpolate(A, B, C)?

Thank you very much for your help

Answer 2

[Sulaymon Eshkabilov]

There are several fucntions of MATLAB which can be used in such exercises. They are interp1(), interp2, interpn,pchip, spline, griddedInterpolant, scatteredInterpolant, etc. Here is one example how it can be attained.

% Given:
F =    [0 1  3  1  0  -1 -3.5  -2 -1 0];             % Given signal
T =    [0 1  3  4  7   8  10  11  12 13];            % Time intervals corresponding
% Note that T is non-uniform
% This how to perform interpolation:

dt = 0.2;     % Sampling time interval, s
fs = 1/dt;    % Sampling frequency, Hz
ti = T(1):1/fs:T(end);  % Time 
FF = griddedInterpolant(T, F);    % Grid interpolation w.r.t ti's time intervals
fi = FF(ti);
figure(1)
plot(T,F, 'ro-', 'DisplayName', 'Data')
hold on
plot(ti, fi, 'k-.', 'Linewidth',2,'DisplayName', 'Interpolated data')
grid on
legend show
hold off

%% Uniform time interval: Note the difference

% Given:
F =    [0 1  3  1  0  -1 -3.5  -2 -1 0]; % Given signal
T =    0:9;                              % Time intervals are uniform
% Note that T can be also non-uniform
% This how to perform interpolation:

dt = 0.2;    % Sampling time interval, s
fs = 1/dt;    % Sampling frequency, Hz
ti = T(1):1/fs:T(end);  % Time 
% FF = griddedInterpolant(T, F);    % Grid interpolation w.r.t ti's time intervals
% fi = FF(ti);
fi = interp1(T,F, ti, 'makima');
figure(2)
plot(T,F, 'ro-', 'markerfacecolor', 'c', 'DisplayName', 'Data')
hold on
plot(ti, fi, 'k-.', 'linewidth',2, 'DisplayName', 'Interpolated data')
grid on
hold off
legend show

[Valeriy]

@Sulaymon Eshkabilov, thanks a lot for your detailed comment. But, as it follows from second part of my previous comment, so I'm repeating:

Second, sampling of my signal realized with constant, uniform spatial frequency, so how to understand interpolation of the original signal to uniform grid? It is already in uniform grid.

Is shape function from this publication shape of my periodic signal?

Probably, to clarify my question, I have to recover shape of the function, which I register as following:

All data sampled with constant, regular, known sampling rate. But how to process them according reference you gave me?

Thank you very much for your help

Answer 3

[Sulaymon Eshkabilov]

As understood your question correctly, is that you are trying to get the shape function which is the envelope thata can be computed by the following fcn of MATLAB - envelope() - see this example:

% DATA: your data 
fm  = 5;          % Freq  1
fc = 100;         % Freq 2
dt = 0.001;       % Sampling time interval, s
fs = 1/dt;        % Sampling frequency, Hz
t= 0:1/fs:1;      % Time 
Am=3;
Ac = 3;
Ka = .75;
m = Am*cos (2*pi*fm*t); 
S = Ac*(1 + Ka*m).* cos(2*pi*fc*t);
%%

% ENvelope is computed here:
[EUp, ELow]=envelope(S);   % Upper and Lower envelopes are computed here
plot(t,S)
hold on
plot(t, EUp, 'r-'),
plot(t, ELow, 'r-')
legend('Data', 'Envelope')
grid on
xlabel('Time, [s]')
ylabel('Signal')

[Valeriy]

@Sulaymon Eshkabilov, thank you, what I'm doing now is something similar to "envelope" you are proposing. I'm calculating it, using

movmax(Signal,nAvrg);
movmin(Signal,nAvrg);

but for cases as last my illustration, nAvrg is huge, several thousends samples. So I significantly lost resolution and this is the reason of my question here and attempts to apply results of the article you are referenced. I hope to find signal shape, using smaller amount of the samples.

By the way,

[EUp, ELow]=envelope(S);

is located in what toolbox? My search shows no results.

[Sulaymon Eshkabilov]

The proposed answer is acceptable :) ?

The envelope() fcn is in the signal processing toolbox. This is the fcn, i.e.,:

function [upperEnv,lowerEnv] = envelope(x, n, method)
%ENVELOPE Envelope detector.
%   [YUPPER,YLOWER] = ENVELOPE(X) returns the upper and lower envelopes of
%   the input sequence, X, using the magnitude of its analytic signal.
%   
%   The function initially removes the mean of X and restores it after
%   computing the envelopes.  If X is a matrix, ENVELOPE operates
%   independently over each column of X.
%
%   [YUPPER,YLOWER] = ENVELOPE(X,N) uses an N-tap Hilbert filter to compute
%   the upper envelope of X.
%
%   [YUPPER,YLOWER] = ENVELOPE(X,N,ENVTYPE) specifies the type of envelope
%   to return. The default is 'analytic':
%      'analytic' - returns the analytic envelope via an N-tap FIR filter
%      'rms'      - returns the RMS envelope of X over a sliding window
%                   of N samples.  
%      'peak'     - returns the peak envelope of the signal using a spline
%                   over local maxima separated by at least N points.
%
%   ENVELOPE(...) without output arguments plots the signal and the upper
%   and lower envelope.
%
%   % Example 1:
%   %   Plot the analytic envelope of a decaying sinusoid.
%   x = 1 + cos(2*pi*(0:999)'/20).*exp(-0.004*(0:999)');
%   envelope(x);
%
%   % Example 2:
%   %   Plot the analytic envelope of a decaying sinusoid using a filter
%   %   with 50 taps.
%   x = 1 + cos(2*pi*(0:999)'/20).*exp(-0.004*(0:999)');
%   envelope(x,50);
%
%   % Example 3:
%   %   Compute the moving RMS envelope of an audio recording of a train
%   %   whistle over every 150 samples.
%   load('train');
%   envelope(y,150,'rms');
%   
%   % Example 4:
%   %   Plot the upper and lower peak envelopes of a speech signal
%   %   smoothed over 30 sample intervals.
%   load('mtlb');
%   envelope(mtlb,30,'peak');
%   
%   See also: HILBERT RMS MAX MIN.

%   Copyright 2015-2018 The MathWorks, Inc.

%   References:
%     [1] Alan V. Oppenheim and Ronald W. Schafer, Discrete-Time
%     Signal Processing, 2nd ed., Prentice-Hall, Upper Saddle River, 
%     New Jersey, 1998.

narginchk(1,3);
nargoutchk(0,2);

if nargin > 2
    method = convertStringsToChars(method);
end

validateattributes(x,{'single','double'},{'2d','real','finite'});

needsTranspose = isrow(x);
if needsTranspose
  x = x(:);
end

if nargin>1
  validateattributes(n,{'numeric'},{'integer','scalar','positive'}, ...
    'envelope','N',2);
end

% allow for 'peaks' (common typo) although it should be 'peak' envelope
if nargin>2
  method = validatestring(method,{'analytic','rms','peaks'});
else
  method = 'analytic';
end

if strcmpi(method,'peaks')
  % no need to remove DC bias from peak finding algorithm
  [yupper, ylower] = envPeak(x,n);
else
  % remove DC offset
  xmean = mean(x);
  xcentered = bsxfun(@minus,x,xmean);
  
  % compute envelope amplitude
  if nargin==1
    xampl = abs(hilbert(xcentered));
  elseif strcmpi(method,'analytic')
    xampl = envFIR(xcentered,n);
  elseif strcmpi(method,'rms')
    xampl = envRMS(xcentered,n);
  end
  
  % restore offset
  yupper = bsxfun(@plus,xmean,xampl);
  ylower = bsxfun(@minus,xmean,xampl);
end

if nargout==0
  plotEnvelope(x,yupper,ylower,method);
elseif needsTranspose
  upperEnv = yupper';
  lowerEnv = ylower';
else
  upperEnv = yupper;
  lowerEnv = ylower;
end

function plotEnvelope(x,yupper,ylower,method)
% plot each signal with a muted default color
% plot each envelope with the default color
colors = get(0,'DefaultAxesColorOrder');
for i=1:size(x,2)
  lineColor = colors(1+mod(i-1,size(colors,1)),:);
  if i==1
    hLine = plot(x(:,1));
    hAxes = ancestor(hLine,'axes');
    axesColor = get(hAxes,'Color');
    set(hLine,'Color',(axesColor+lineColor)/2);
  else
    line(1:size(x,1),x(:,i),'Color',(lineColor+axesColor)/2);
  end
  line(1:size(x,1),yupper(:,i),'Color',lineColor);
  line(1:size(x,1),ylower(:,i),'Color',lineColor);
end

% add legend only when one signal is present
if size(x,2)==1
  legend(getString(message('signal:envelope:Signal')), ...
         getString(message('signal:envelope:Envelope')));
end

% lookup and plot title string
catStrs = {'Analytic','RMS','Peak'};
catStr = catStrs{strcmp(method,{'analytic','rms','peaks'})};
titleStr = getString(message(['signal:envelope:' catStr 'Envelope']));
title(titleStr);


function y = envFIR(x,n)

% construct ideal hilbert filter truncated to desired length
fc = 1;
t = fc/2 * ((1-n)/2:(n-1)/2)';

hfilt = sinc(t) .* exp(1i*pi*t);

% multiply ideal filter with tapered window
beta = 8;
firFilter = hfilt .* kaiser(n,beta);
firFilter = firFilter / sum(real(firFilter));

% apply filter and take the magnitude
y = zeros(size(x),'like',x);
for chan=1:size(x,2)
  y(:,chan) = abs(conv(x(:,chan),firFilter,'same'));
end

% compute RMS 
function y = envRMS(x,n)
y = signal.internal.builtin.movrms(x,n,'same');

function [yupper,ylower] = envPeak(x,n)

% pre-allocate space for results
nx = size(x,1);
yupper = zeros(size(x),'like',x);
ylower = zeros(size(x),'like',x);

% handle default case where not enough input is given
if nx < 2
  yupper = x;
  ylower = x;
  return
end

% compute upper envelope
for chan=1:size(x,2)
  if nx > n+1
    % find local maxima separated by at least N samples
    [~,iPk] = findpeaks(double(x(:,chan)),'MinPeakDistance',n);
  else
    iPk = [];
  end
  
  if numel(iPk)<2
    % include the first and last points
    iLocs = [1; iPk; nx];
  else
    iLocs = iPk;
  end

  % smoothly connect the maxima via a spline.
  yupper(:,chan) = interp1(iLocs,x(iLocs,chan),(1:nx)','spline');
end

% compute lower envelope
for chan=1:size(x,2)
  if nx > n+1
    % find local minima separated by at least N samples
    [~,iPk] = findpeaks(double(-x(:,chan)),'MinPeakDistance',n);
  else
    iPk = [];
  end
  
  if numel(iPk)<2
    % include the first and last points
    iLocs = [1; iPk; nx];
  else
    iLocs = iPk;
  end
  
  % smoothly connect the minima via a spline.
  ylower(:,chan) = interp1(iLocs,x(iLocs,chan),(1:nx)','spline');
end

[Valeriy]

@Sulaymon Eshkabilov, I appreciate a lot your help, thank you.

The proposed answer is acceptable :) ?

Your first answer with the link to an article is most close to being acceptable. However, one crucial point remains: I can't realize their proposal due to the lack of necessary mathematical level. So this is the reason why I'm asking for your help.

The method of how to answer my question by envelope (or equivalent) function I have realized before asking this question.

Thank you for your explanation about transferring my data to a uniform grid, but my data are uniformly spaced. The next stage of Algorithm 1 from the cited article is 2. Compute the Fourier coefficients of f in theta space. In my case, a signal changes amplitude along many sampling periods. How long should signal data to be taken if its amplitude is changed over thousands of periods? In fact, my question has to be reformulated, how big has to be N from the first point of Algorithm 1?

Is Matlab realization of the singular value decomposition (SVD) necessary in the fifth stage of Algorithm 1?

Answer 4

[Valeriy]

@Sulaymon Eshkabilov, I understand that answer on my question, using methods of cited article (https://royalsocietypublishing.org/doi/10.1098/rsta.2015.0194) is not too simple and perhaps it is behind your scope. But I'm sure that there are some peoples, members of this community with background, which well corresponds to knowledge, which is necessary to decipher proposed methods.

Do you know the way, methods, how to ask them, how to pass this question to correct audience? Make extended copy of Algorithm 1 from article?

Thanks for your help and proposals