reproduces the audio signal after sampling!!!

One of the earliest extensions of this theorem was stated by Shannon himself in his 1949 paper, which says that if x(t)and its first (M - 1) derivatives are available, then uniformly spaced samples of these, taken at the reduced rate of M times, are sufficient to reconstruct x(t) . This result will be referred to as the derivative sampling theorem in this paper. I'm working with M=3 and having trouble designing a synthesis filter with an impulse response of the following form: sinc(t)^3; t.sinc(t)^3; t^2.sinc(t)^3. please help me?

[Paul]

Hi Vu Hoang Son,

Is the goal to define three filters in the frequency domain, one that has impulse response sinc(t)^3, one that has impulse response t*sinc(t)^3, and one that has impulse response t^2*sinc(t)^3?

The ; is confusing me.

[Vu Hoang Son]

You're right! My problem is that I need to design three synthesis filters to recover the signal based on sampling the signal and sampling the first and second derivatives of the signal. Can you give me a solution?

Answer 1

[Paul]

Hi Vu,

I wasn't familiar with derivative sampling, or generalized sampling, until I saw your post. After a bit of reading, maybe I'm still not because I can't get it to work for N = 3, though it does for N = 1 and N = 2.

Define a signal for analysis, where it's understood that at t = 0 x(t) will be defined by its limit. Note that x(t) is continuous.

syms w t real
x(t) = (cos(sym(pi)*t)-1)^2/t^4/sym(pi)^4

x(t) = 

figure
fplot(real(x(t)),[-10 10]),xlabel('t');ylabel('x(t)')

Get the Fourier transform of x(t)

X(w) = fourier(rewrite(x(t),'exp'),t,w)

X(w) = 

Convert to Hz

syms f real
X(f) = X(2*sym(pi)*f);
imag(X(f))

ans = 

0

figure
fplot(real(X(f))),xlabel('f'),ylabel('X(f)')

Show that X(f) is bandlimited, i.e., X(f) = 0 for abs(f) > 1.

syms f0 real
assume(f0 > 1);
simplify(X(f0))

ans = 

0

assume(f0 < -1);
simplify(X(f0))

ans = 

0

First derivatve of x(t)

dx(t) = simplify(diff(x(t),t))

dx(t) = 

figure
fplot(dx(t),[-10 10]);xlabel('t'),ylabel('dx/dt')

Second derivative of x(t)

d2x(t) = simplify(diff(dx(t),t))

d2x(t) = 

figure
fplot(d2x(t),[-10 10]),xlabel('t'),ylabel('d^2x/dt^2')

Get the values of all three at t = 0

limit([x(t), dx(t), d2x(t)],t,0)

ans = 

Create functions for numerical evaluation all three.

temp = matlabFunction((x(t)));
xfunc = @(t) fillmissing(temp(t),'constant',1/4);

temp = matlabFunction((dx(t)));
xdfunc = @(t) fillmissing(temp(t),'constant',0);

temp = matlabFunction((d2x(t)));
xd2func = @(t) fillmissing(temp(t),'constant',-(pi^2)/12);

Define the bandwidth, the Nyquist frequency, and the Nyquist sampling period:

fband = 1;
fnyquist = 2*fband;
delta = 1/fnyquist;

Derivative interpolation for N = 1,2, and 3

for N = 1:3

    Ts = delta*N;
    % based on plots above, it looks like sampling out to 10 seconds should
    % be sufficient.
    nmax = round(10/Ts);
    n = -nmax:nmax;
    
    % sample the signal and its derivatives
    xsamp   = xfunc  (n*Ts);
    xdsamp  = xdfunc (n*Ts);
    xd2samp = xd2func(n*Ts);
    
    % reconstruction
    xval = cat(3,xsamp,xdsamp,xd2samp);
    
    k = reshape(0:N-1,1,1,[]);
    
    tval = (-10:.01:10).';
    
    M = ((tval - N*delta*n).^k).*(sinc((tval - N*delta*n)/(N*delta)).^N)./factorial(k);
    xr = sum(M.*xval(:,:,1:N),[2 3]);
    
    figure
    plot(tval,xfunc(tval))
    hold on
    plot(n*Ts,xsamp,'o')
    plot(tval,xr)
    legend('x(t)','x(n*Ts)','xr(t)')
    xlim([-10 10])
    xlabel('t'),ylabel('x')

end

I'm at a loss as to why it doesn't work for N = 3, though the reconstructed signal does go through the samples. Don't know if it's a coding error, or maybe x(t) doesn't satisfy some requirement for this type of sampling, or something else.

Maybe someone else (@David Goodmanson, @Matt J) will have some insight.

[Vu Hoang Son]

Hi @Paul ! Thank you for your interest in my question. Currently I don't know the exact answer. If you have any ideas to solve this problem please help me. Sincerely thank

[Paul]

At present, I'm still stumped why the code seems to give the expected result for N = 1 and N = 2, but not N = 3 ...

[Paul]

The code doesn't give the expected result for N = 3 because the reconstruction equation shown in this comment is incorrect. The equation in that comment is the same as equation (1) in this paper by Linden and Abramson. However, the authors subsequently posted a correction to that paper as was identified by a user at this link at dsp.stackexchange.com .

For the case at hand with N = 3 (which is R = 2 in the notation of Linden and Abramson), the third term in the summation does indeed include the additional term (z) derived by @David Goodmanson in this answer.

[Paul]

Here is the example with code based on the Linden and Abramson correction, and using their notation.

Define the example signal and its Fourier transform

syms w t real
x(t) = (cos(sym(pi)*t)-1)^2/t^4/sym(pi)^4;  

X(w) = fourier(rewrite(x(t),'exp'),t,w);
syms f real
X(f) = simplify(X(2*sym(pi)*f));

Verify that x(t) is bandlimited with X(f) = 0 for f > 1

W = 1;

syms f0 real
assume(f0 > W);
simplify(X(f0))

ans = 

0

assume(f0 < -W);
simplify(X(f0))

ans = 

0

Define the maximum value of R to explore. R+1 is the scale factor on the Nyquist sampling period

Rmax = 6;

Define anonymous functions to evaluate x(t) and its derivatives

for r = 0:Rmax
    dx = simplify(diff(x(t),t,r));
    lim = double(limit(dx,t,0));
    temp = matlabFunction(dx);
    xf{r+1} = @(t) fillmissing(temp(t),'constant',lim);
end

Values at which to reconstruct x(t)

tval = (-10:.01:10).';

Reconstruction for different sampling factors

for R = 0:Rmax

    xr = zeros(size(tval));

    % sampling period
    tau = (R+1)/(2*W);

    % sampling out to 10 seconds should be sufficient.
    kmax = round(10/tau); % k
    k    = -kmax:kmax;    % k

    sincval = sinc((tval - k*tau)/tau).^(R+1);
    
    % Incrementally sum to form the reconstruction xr(t)
    for r = 0:R

        % compute eta based on samples of x(t) and its derivatives
        eta = zeros(size(k));
        jj = r;
        for ii = 0:jj
            eta = eta + nchoosek(jj,ii)*(pi/tau)^(jj-ii)*gam(R+1,jj-ii,t)*xf{ii+1}(k*tau);
        end

        xr = xr + sum((tval - k*tau).^r.*sincval/factorial(r).*eta,2);

    end
    
    % Plot the true x(t) overlayed with its samples and the reconstruction xr(t)
    figure
    plot(tval,xf{1}(tval))
    hold on
    plot(k*tau,xf{1}(k*tau),'o')
    plot(tval,xr)
    legend('x(t)','x(n*Ts)','xr(t)')
    xlim([-10 10])
    xlabel('t'),ylabel('x')
    title("R = " + string(R))

end

function g = gam(alpha,beta,t)

if mod(beta,2) == 0  % beta even

    g = double(limit(diff((t/sin(t))^alpha,t,beta),t,0));

else

    g = 0;

end

end

[Vu Hoang Son]

@Paul Are you interested in signal processing? Can you help me apply it to voice signals?

[Paul]

Yes, I have a hobby-like interest in signal processing. I'd only respond to Matlab questions posted here on Answers.

[Vu Hoang Son]

To apply to real voice signals is not simple. I'm completely stuck....:((

Answer 2

[David Goodmanson]

Hello Vu & Paul,

< I took a brief look at the Vaidyanathan.paper and although the decimation process seems clear enough, the 'interpolation' and recombination might differ from straight addition of the three components, in which case the conclusion here about the need for an exra term might not apply. >

I believe I have the result corresponding to the expression containing sinc functions shown in Vu's last comment. (the OP expression). And I believe that for N = 3 and larger, that expression is incorrect, hence Paul's result. I use cap Delta --> ts (sampling time) , N*delta --> Nts and f(t) --> y(t) below.

All sums over n run from -inf to inf.

The expression is, for known sample values of y,y',y''.

y(t) = sum{n}
[ y(n*Nts) + y'(n*Nts)*(t-n*Nts) + ( z + y''(n*Nts) )*(t-n*Nts)^2/2 ]
                                                               *sinc((t-n*Nts)/Nts)^3

but z represents a missing term which is

(pi/Nts)^2*y(n*Nts)

and involves known values of y, not of y''. Once that is included the results look fairly good.

To simplify things, for both the signal and all sinc functions I used functions all of whose fourier transfoms are in the same interval, -f0<=f<=f0, which I will call condition F. Using different frequency intervals for the signal and for the sinc functions probably works in some circumstances but could be an additional complication.

Hand waving here, but If

transform of     t^p*sinc(a*t)     is in       -f0<=f<=f0       condition F 
then
transform of    (t^p*sinc(a*t))^q    is in     -q*f0<=f<=q*f0 
and
transform of    (t^p*sinc(a*t/q))^q  is in       -f0<=f<=f0        back to condition F

So the arguments of sinc() are adjusted accordingly to compensate. (This happens automatically in the OP expression by having Nts in the denominator of the sinc argument).

I didn't really doubt Paul's results with symbolics but I wanted to do things differently for comparison so the calculation is strictly numeric, including calculation of y' and y'' for the input signal.

If lots of terms are used in the summation (nmax below) then eventually the sampling points extend beyond the time window. When that happens the results go bad, probably due to an aliasing effect that I have not tried to figure out. The code keeps the o's within the window.

A derivation of the z term is shown below the plots, with code after that. I don't want to write an essay here, but since people have been referring to published results I'll take up more space than usual.

To start, for N=1 let ts be the sampling interval, and let u0(t) be a continuous function whose transform satisfies condition F above and for which

u0(0) = 1     u0(n*ts) = 0  n ~=0

or which is the same thing, u0 has the orthogonality property

u0((m-n)*ts) = delta(m,n)                                          (1)

Create a continuous function using known values at the sampling points

y(t) = sum{n} y(n*ts)*u0(t-n*ts)     

Evaluate this @ t = m*ts, and with (1) show that y(m*ts) = y(m*ts), equal on both sides, which proves the sampling theorem. Here

u0(t) = sin(pi*t/ts)/(pi*t/ts) = sinc(t/ts)

works, with the Nyquist condition

2*f0*ts = 1.

(I wish that the nomenclature sinc(x) = sin(pi*x)/(pi*x) had not taken over because the mathematical definition sinc(x) = sin(x)/x predates the signal processing definition, but that's the way it goes).

Now for N = 3, denote

N*ts = Nts

and suppose there are three functions (shown on the right) that all meet condition F and for which

  u0((m-n)*Nts) = delta(m,n)                      u0(t) =           sinc(t/Nts)^3
 u1'((m-n)*Nts) = delta(m,n)                      u1(t) =         t*sinc(t/Nts)^3
u2''((m-n)*Nts) = delta(m,n)                      u2(t) =   (t^2/2)*sinc(t/Nts)^3

As mentioned before, in order to to stay with condition F, sinc(t/Nts)^3 undersamples the signal by a factor N=3. But this is made up by the conditions on y and y'.

for u0 , evaluating @ t = m*Nts similarly to what was done before,

y(t) = sum{n} y(n*Nts)*u0(t-n*Nts)           -->  y(m*Nts) = y(m*Nts)
                 since u0 meets the orthogonality property  

For the others,take derivatives, evaluate @ t = m*Nts

  y(t) = sum{n}  y'(n*Nts)* u1(t-n*Nts)
 y'(t) = sum{n}  y'(n*Nts)*u1'(t-n*Nts)    -->  y'(m*Nts) = y'(m*Nts)
                 since u1' meets the orthogonality property  

  y(t) = sum{n} y''(n*Nts)*  u2(t-n*Nts)
y''(t) = sum{n} y''(n*Nts)*u2''(t-n*Nts)    -->  y''(m*Nts) = y''(m*Nts)
                    since u2'' meets the orthogonality property

These all match up, so the function that is the sum of these three,

y(t) = sum{n} [ (y(n*Nts)*u0(t-n*Nts) + y'(n*Nts)*u1(t-n*Nts) + y''(n*Nts)*u2(t-n*Nts) ]
                                                                                   (4)

as in the OP expression, works. Except that it doesn't. The problem is that you have to look at y and its derivatives etc. after the sum of the three is done. For example,

y''(t) = sum{n} [(y(n*Nts)*u0''(t-n*Nts)+ y'(n*Nts)*u1''(t-n*Nts) + y''(n*Nts)*u2''(t-n*Nts)]
                                                                                   (5)

and this will only work if u0'' = u1'' = 0 for all sample points. The same considerations apply to the y(t) and y'(t) equations, and overall what is needed is for all the 'nondiagonal' contributions

u0' = u0'' = u1 = u1'' = u2 = u2' = 0

for all sample points. As it turns out all of these are 0 except for u0'' which as you can check is

u0''(0) = -A     u0''(n*Nts) = 0  n ~=0           where A has the value  (pi/Nts)^2

which is the same thing as the orthogonality relation

u0''((m-n)*Nts) = -A*delta(m,n)             

Therefore evaluating (5) @ t = m*Nts gives an extra unwanted term:

y''(m*Nts) = -A*y(m*Tns) + y''(m*Nts)    ?

However this can be canceled out by introducing an extra term into (4) of

z*u2(t-n*Nts) = [A*y(n*Tns)]*u2(t-n*Nts) 

(which does not interfere with anything else since u2 = u2' = 0 at at all sampling points). This is the z term which gives the correct result for y''. Since

u2(t-n*Nts) = ((t-n*Nts)^2/2)*sinc((t-n*Nts)/Nts)^3,

this jibes with the rewrite of the OP expression at the top of this answer.

precis = [];                      % check precision of the solution
t = -30:.001:30;
f0 = 1;
y = sinc(2*f0*t/4).^4;            % signal
y1 = splinder(t,y);               % 1st deriv
y2 = splinder(t,y1);              % 2nd deriv
% N = 1;
ts = 1/(2*f0);                    % nyquist
N = 1
Nts = N*ts;
nmax = floor(24/Nts);
n = (-nmax:nmax)';
yn0 = spline(t,y,n*Nts);         % sample points, row vec
[tmat nNmat] = meshgrid(t,n*Nts);
u0 = sinc((tmat-nNmat)/Nts);
yr = sum(yn0.*u0);               % reconstructed
figure(14); subplot(2,1,1)
plot(n*Nts,yn0,'o',t,y,t,yr)
grid on
legend('sample points', 'y','yr')
title(['N  =  ' num2str(N)])
precis = [precis; max(abs(y-yr))];
% N = 2;
ts = 1/(2*f0);                     % nyquist
N = 2
Nts = N*ts;
nmax = floor(24/Nts);
n = (-nmax:nmax)';
yn0 = spline(t,y,n*Nts);          % sample points, row vec
yn1 = spline(t,y1,n*Nts);         % 1st deriv sample points
[tmat nNmat] = meshgrid(t,n*Nts);
u0 = sinc((tmat-nNmat)/Nts).^2;
u1 = (tmat-nNmat).*sinc((tmat-nNmat)/Nts).^2;
yr = sum(yn0.*u0 + yn1.*u1);
figure(14); subplot(2,1,2)
plot(n*Nts,yn0,'o',t,y,t,yr)
grid on
legend('sample points', 'y','yr')
title(['N  =  ' num2str(N)])
precis = [precis; max(abs(y-yr))];
% N = 3;
ts = 1/(2*f0);                     % nyquist
N = 3
Nts = N*ts;
nmax = floor(24/Nts);
n = (-nmax:nmax)';
yn0 = spline(t,y,n*Nts);          % sample points, row vec
yn1 = spline(t,y1,n*Nts);         % 1st deriv sample points
yn2 = spline(t,y2,n*Nts);         % 2ndt deriv sample points
[tmat nNmat] = meshgrid(t,n*Nts);
u0 = sinc((tmat-nNmat)/Nts).^3;
u1 = (tmat-nNmat).*sinc((tmat-nNmat)/Nts).^3;
u2 = ((tmat-nNmat).^2/2).*sinc((tmat-nNmat)/Nts).^3;
yrwith =    sum(yn0.*u0 + yn1.*u1 +(yn2 + (pi/Nts)^2*yn0).*u2);
yrwithout = sum(yn0.*u0 + yn1.*u1 +(yn2 +      0        ).*u2);
figure(16); subplot(2,1,1)
plot(n*Nts,yn0,'o',t,y,t,yrwith)
grid on
legend('sample points', 'y','yr')
title(['N  =  ' num2str(N) '  wih ''z'' term'])
figure(16); subplot(2,1,2)
plot(n*Nts,yn0,'o',t,y,t,yrwithout)
grid on
legend('sample points', 'y','yr')
title(['N  =  ' num2str(N) '  without ''z'' term'])
precis = [precis; max(abs(y-yr))]  % precsion of the result for N = 1,2,3

function dydx = splinder(x,y,x1)

% Differentiation of y(x) by cubic spline 
% array x must be sorted in ascending order
% splinder(x,y,x1) is the derivative of y(x) at the points in vector x1.
% splinder(x,y) sets x1 = x and is faster than splinder(x,y,x).
% Output format (row or column vector) matches x1.
%
% dydx = splinder(x,y,x1)

% David Goodmanson
% function changed to this name july 2015

if nargin == 2;  x1 = x;  end;
[brk,c] = unmkpp(spline(x,y));
% c has cubic coeffs in col 1, const coeffs in col 4
[nrow ncol] = size(c);
for k = 1:ncol
  c(:,k) = (ncol-k)*c(:,k);
end
c(:,end) = [];

if nargin == 2 & length(x) >=4    % quick method when x1 = x.
  lastone = polyval(c(end,:),brk(end)-brk(end-1));
  dydx = [c(:,end); lastone] ;
  if size(x,1) == 1
    dydx = dydx.';
  end
else                               % x1 different from x
  dydx = ppval(mkpp(brk,c),x1);
end

[Paul]

Hi David,

Thanks for posting this. Do you have a reference that shows the additional term? I'll wait for your "more to follow."

As this topic was, and still is, new to me, I've been trying to find basic, text-book-like material with an example or two, but have not been successful, at least with the time and effort that I've been willing to invest.

I did find this paper by Linden and Abramson that shows the same equation posted by @Vu Hoang Son (see eq (1)). Hence, I've been concerned that my function doesn't satisify some assumptions that lead to their result, though I've been unable to convince myself that such a concern is actually warranted.

I've also been wondering if the N = 3 case converges really slowly and many, many more samples are needed, small as they may be. I've played around with increasing nmax, but that didn't really seem to help.

This other paper by Jagerman and Fogel is a more math-y and has some additional, interesting Theorems. But, unfortunately, it only shows the reconstruction expression using samples up to the first derivative, whereas we want to also use samples of the second derivative.

Another paper by Fogel shows the general case, but, unfortunately for me, doesn't show the closed form expression for the particular case of using samples of the 0th, first, and second derivatives.

[David Goodmanson]

Hi Paul, I can get back to this in a couple of hours (I worked on this intently and realized I'm starving; I'm sure you know how that can go when you're caught up with equations) but I derived the missing term and am convinced that for N=3 the problem is not with slow convergence, that too many tterms in the sum eventually foul things up as I metioned and that the problem is due to the missing term.

[Vu Hoang Son]

Hi David & Paul,

Thank you for your interest in my topic. @David Goodmanson please explain to me why we need to add z when N>3? I have read some publications where they did derivative sampling with N=2 and then generalized it to N=3 and larger, especially with fractional derivative sampling. Some of that work is based on filter banks (e.g., "Filterbank Reconstruction of Bandlimited Signals from Nonuniform and Generalized Samples" by Yonina C. Eldar and Alan V. Oppenheim; "On the reconstruction of derivative sampling method of band-limited signal" by Chien-Cheng Tseng, Su-Ling Lee. Or some works based on the fractional Furie transform (FRFT) domain (e.g. "Second Order Derivative Sampling and Its Application into Image Scaling" by Yong Li; "Higher-Order Derivative Sampling Associated with Fractional Fourier Transform" by Rui-Meng Jing, Qiang Feng, Bing- Zhao Li ). I read it but really didn't understand much and couldn't apply it to my problem, which is speech signal processing with 3rd sampling and 4th derivative. It took me a long time but no progress yet! @Paul @David Goodmanson please help me. If you have any ideas please suggest me!!!!

[David Goodmanson]

Hello V&P, a modified answer covers the derivation of the extra z term.

[Vu Hoang Son]

Hi David, Please help me consider the input signal as a voice signal with a frequency range of 0.3-3.4kHz. What will the results be?

[Vu Hoang Son]

From the mathematical expression shown, I want to model it in the system as shown. Can you help me @David Goodmanson@Paul?

[Paul]

@Vu Hoang Son

Do you have a reference for the expression you showed in this comment?

[Vu Hoang Son]

@Paul I have referred to the document "Classical Sampling Theorems in the Context of Multirate and Polyphase Digital Filter Bank Structures" by P.P. Vaidyanathan. The exact expression above I took in a document in Russian. If needed, give me your email address and I will send it to you

[Vu Hoang Son]

@David Goodmanson@PaulSee if this document is helpful. I read it but it's too difficult

[Paul]

@David Goodmanson

Are there any conditions that can be imposed on y(t) such that the additional term, z, is not needed? Presumably such a condition would actually be imposed on d2y(t)/dt2.

[David Goodmanson]

If you look at backpedal note I inserted at the top of the answer, I don't know enough signal processing to know how things are recombined, but if it not just straight addition then maybe the whole issue is moot.

Answer 3

[Paul]

The Fourier transform of each of those impulse responses can be found as follows. See also this thread for related discussion.

syms t w real
f(t) = sinc(t)^3

f(t) = 

Matlab can't find the Fourier transform

F(w) = fourier(f(t),t,w)

F(w) = 

But, it can if we rewrite f(t) in terms of exp

f(t) = rewrite(f(t),'exp')

f(t) = 

F(w) = simplify(fourier(f(t),t,w))

F(w) = 

figure
fplot(F(w),[-20 20])

Then, for t*f(t)

F1(w) = simplify(fourier(t*f(t),t,w))

F1(w) = 

real(F1(w))

ans = 

0

figure
fplot(imag(F1(w)),[-20 20])
ylim([-0.2 0.2])

Similarly

F2(w) = simplify(fourier(t^2*f(t),t,w))

F2(w) = 

figure
fplot(F2(w),[-20 20])
ylim([-0.06 0.06])

[Vu Hoang Son]

Thank you very much. My problem is applying these 3 filters in speech signal processing. The input is a speech signal, which is then sampled including the signal, the first derivative of the signal, and the second derivative of the signal. Can using the above 3 filters restore the original signal? How to simulate it in matlab? Help me please!

[Vu Hoang Son]

This is the formula that represents the problem I am talking about. k-th derivative of the signal. - Sampling period according to Nyquist-Shannon sampling theorem