# reproduces the audio signal after sampling!!!

124 views (last 30 days)
Vu Hoang Son on 17 Apr 2024
Commented: Vu Hoang Son on 8 Jun 2024 at 21:08
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 on 19 Apr 2024
Edited: Paul on 19 Apr 2024
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 on 20 Apr 2024
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?

Paul on 26 Apr 2024
Edited: Paul on 27 Apr 2024
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.
Paul on 8 Jun 2024 at 20:35
Yes, I have a hobby-like interest in signal processing. I'd only respond to Matlab questions posted here on Answers.
Vu Hoang Son on 8 Jun 2024 at 21:08
To apply to real voice signals is not simple. I'm completely stuck....:((

David Goodmanson on 25 May 2024 at 9:38
Edited: David Goodmanson on 29 May 2024 at 7:54
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 on 30 May 2024 at 3:16
Edited: Paul on 30 May 2024 at 11:47
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 on 30 May 2024 at 15:49
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.

Paul on 20 Apr 2024
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 on 21 Apr 2024
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 on 21 Apr 2024
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

### Categories

Find more on Multirate Signal Processing in Help Center and File Exchange

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!