reproduces the audio signal after sampling!!!

Question

Vu Hoang Son 2024-4-17

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2108126-reproduces-the-audio-signal-after-sampling

评论： Vu Hoang Son 2024-6-8

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?

2 个评论
显示无隐藏无

Paul 2024-4-19

编辑：Paul 2024-4-19

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 2024-4-20

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 2024-4-26

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2108126-reproduces-the-audio-signal-after-sampling#answer_1448521

编辑：Paul 2024-4-27

在 MATLAB Online 中打开

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.

7 个评论
显示 5更早的评论隐藏 5更早的评论

Paul 2024-6-8

编辑：Paul 2024-6-8

在 MATLAB Online 中打开

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

Paul 2024-6-8

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

Vu Hoang Son 2024-6-8

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

请先登录，再进行评论。

Answer 2

David Goodmanson 2024-5-25

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2108126-reproduces-the-audio-signal-after-sampling#answer_1463276

编辑：David Goodmanson 2024-5-29

在 MATLAB Online 中打开

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

11 个评论
显示 9更早的评论隐藏 9更早的评论

Paul 2024-5-25

编辑：Paul 2024-5-25

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.

Vu Hoang Son 2024-5-26

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

Paul 2024-5-30

编辑：Paul 2024-5-30

@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 2024-5-30

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 2024-4-20

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2108126-reproduces-the-audio-signal-after-sampling#answer_1444901

在 MATLAB Online 中打开

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

2 个评论
显示无隐藏无

Vu Hoang Son 2024-4-21

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 2024-4-21