x(t) =

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.

25 views (last 30 days)

Show older comments

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

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)

Convert to Hz

syms f real

X(f) = X(2*sym(pi)*f);

imag(X(f))

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))

assume(f0 < -1);

simplify(X(f0))

First derivatve of x(t)

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

figure

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

Second derivative of x(t)

d2x(t) = simplify(diff(dx(t),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)

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.

Paul
on 8 Jun 2024

David Goodmanson
on 25 May 2024

Edited: David Goodmanson
on 29 May 2024

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

David Goodmanson
on 30 May 2024

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

Matlab can't find the Fourier transform

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

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

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

figure

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

Then, for t*f(t)

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

real(F1(w))

figure

fplot(imag(F1(w)),[-20 20])

ylim([-0.2 0.2])

Similarly

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

figure

fplot(F2(w),[-20 20])

ylim([-0.06 0.06])

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

Start Hunting!