Skip to main content ↑ ↓

Notice

Please note that most of the software linked on this forum is likely to be safe to use. If you are unsure, feel free to ask in the relevant topics, or send a private message to an administrator or moderator. To help curb the problems of false positives, or in the event that you do find actual malware, you can contribute through the article linked here.
Topic: "Unique" band limited function (Digital Show and Tell) (Read 6883 times) previous topic - next topic
0 Members and 1 Guest are viewing this topic.

"Unique" band limited function (Digital Show and Tell)

I recently watched monty's excellent video, "Digital Show and Tell", but found one of the claims confusing. He asserts that there is a unique band-limited function that fits the data points. Thus, higher sampling rates add no additional information in the lower frequencies.

This doesn't match my understanding of the Fourier Transform. I learned about the DFT from a background in mathematics---I have never used it for signal analysis.

To make my confusion concrete, run this program in octave:
function f  = interp_fft(y)
  n = size(y,2);
  F = 0:n-1;
  P = fft(y);
  f = @(x) sum(arrayfun(@(f,p) p*exp(2*pi*i*f*x/n), F, P))/n;
endfunction
n=32
x_dot=0:n-1;
x_all=0:0.01:n;
N = rand(1,n);
M = N(1:n/2);
f = interp_fft(N);
g = interp_fft(M);
y = arrayfun(f, x_all);
z = arrayfun(g, x_all);
plot(x_dot,N,'ro',x_all,y,'b-',x_all,z,'k-');

In the resulting graph, we have random data points represented by the 'red dots'. I produced two fitting curves. One curve fits the first half only and the second curve fits all the dots. You can see, that despite agreeing on the first dots, between the dots the functions differ.

Define "16" units on the x-access to be 1 second. Then, as far as I understand this, I turned the 16 data point function into frequency powers at 0Hz, 1Hz, 2Hz, ..., 15Hz. The second curve has power at frequencies at 0Hz, 0.5Hz, 1Hz, 1.5Hz, ...15.5Hz.

Wouldn't both functions still be called "band-limited"? So, I have two functions that pass through the first 16 data points ... and they differ.

What gives?

"Unique" band limited function (Digital Show and Tell)

Reply #1 –
Monty is right. The defect in your program is that FFT cannot be used to produce a bandlimited fit. Bandlimited signals only exist on an infinite support, not on a segment. Your FFT wraps the signal in time around F = 0:n-1

"Unique" band limited function (Digital Show and Tell)

Reply #2 –
Quote from: terpstra on
Wouldn't both functions still be called "band-limited"? So, I have two functions that pass through the first 16 data points ... and they differ.

If both functions were properly bandlimited, you'd see that the difference between them still exists in the first half, but it is decaying as 1/t. Monty's point holds for an infinite set of sampling points, not for 2 finite sets.

"Unique" band limited function (Digital Show and Tell)

Reply #3 –
x(t) can be uniquely reconstructed given the sampling theorem with:

where hs is the sinc function.

Here everything is infinite. Practically you'd for example truncate the sinc filter, decrease its cutoff frequency, multiply it with a window function.. and use that filter to reconstruct x(t).
"I hear it when I see it."

"Unique" band limited function (Digital Show and Tell)

Reply #4 –
Sorry for going off-topic but I amazed by people with this type of mathematical knowledge you master. Where did you all learn these things? Did you attend some seriously advanced education or have you been able to learn this or part of it by yourselves? Once again, sorry for being off-topic.

"Unique" band limited function (Digital Show and Tell)

Reply #5 –
Saying that I need to do the FFT over the whole signal (from start of song to end of song or just from the beginning to end of time) doesn't really answer my concern.

Quote from: Alexey Lukin on
Your FFT wraps the signal in time around F = 0:n-1 If both functions were properly bandlimited, you'd see that the difference between them still exists in the first half, but it is decaying as 1/t

I understand, of course, that using a DFT as I've done it is periodic. Thus my input functions already have the infinite support you desire. However, in the real world of a DAC that produces the resulting analogue function, there must be some sort of windowing going on. The DAC is not an oracle having perfect knowledge of the past and future. Thus, it must do something like I did: select a limited number of samples and fill in the gaps from there.

What is your 't'? Is it the number of samples? If so, then I suppose this answers my question, as it means that a DAC can just use 2^16 samples worth of data in its window and you would reconstruct the signal accurate to the noise floor. What metric is the '1/t' decay measured in? Least squares? Does this convergence theorem have a name?

Quote from: xnor on
x(t) can be uniquely reconstructed given the sampling theorem with... Here everything is infinite. Practically you'd ...

I did more-or-less what you propose in my example code. My reconstructed function comes from the inverse FFT / sampling theorem, and I truncated the data at 32 and 16 respectively.

Ultimately, I guess the subtext of my question is concern about what happens exactly when you go from infinite to discrete. How does the approximation of the true signal converge.

"Unique" band limited function (Digital Show and Tell)

Reply #6 –
Quote from: punkrockdude on
Sorry for going off-topic but I amazed by people with this type of mathematical knowledge you master. Where did you all learn these things? Did you attend some seriously advanced education or have you been able to learn this or part of it by yourselves? Once again, sorry for being off-topic.


I am pretty sure most of us have some sort of university background, and often in something like electrical engineering, where learning that level of mathematics is pretty much a must, even if many of us don't get to use it too much in our work. You could probably learn it on your own, from textbooks and online courseware, but there is a significant risk of getting stuck because of some misunderstanding taking you in the wrong direction unless you get feedback.

"Unique" band limited function (Digital Show and Tell)

Reply #7 –
Quote from: punkrockdude on
Sorry for going off-topic but I amazed by people with this type of mathematical knowledge you master. Where did you all learn these things? Did you attend some seriously advanced education or have you been able to learn this or part of it by yourselves? Once again, sorry for being off-topic.


I don't know about everyone else, but I studied computer engineering, and we had to take a course on this stuff where you basically just sat down and proved all the Fourier and sampling theory stuff.  It sounds really complex, but most of the math isn't that hard.  Its mostly algebra and some introductory calculus, and of course if you're interested in audio, fairly easy to relate to things you understand intuitively.  Compared to a lot of EE or programming classes it was actually pretty easy.

"Unique" band limited function (Digital Show and Tell)

Reply #8 –
Ok, so I went a bit further towards trying to understand this.

First, I improved the performance of my curve fitter so I could test the convergence claim. The first octave program below will plot the two fit functions and their difference:

function Y = fit(y,x,m)
  n = size(y,2);
  F = fft(y);
  G = zeros(1,x);
  G(1:m/n:m) = F;
  Y = ifft(G)*x/n;
endfunction
n=16; 
m=32; 
x=1024;
M=rand(1,m);
N=M(1:n);
YM=fit(M,x,m);
YN=fit(N,x,m);
l=(n-1)*x/m;
E=YM(1:l)-YN(1:l);           
plot([1:x/m:x],M,'ro',[1:x],YM,'b-',[1:x],YN,'k-',[1:l],E,'g-');     

The result is pretty interesting. You can see that the error between the two fits actually looks a lot like a sine wave with the same frequency as the sample rate. The magnitude of the error varies quite a bit between runs, but it generally has the same shape.

The next program tests the claim that more samples = converges to correct curve.
function Y = fit(y,x,m)
  n = size(y,2);
  F = fft(y);
  G = zeros(1,x);
  G(1:m/n:m) = F;
  Y = ifft(G)*x/n;
endfunction
function E = error(n,m,x)
  M = rand(1,m);
  N = M(1:n);
  YM = fit(N,x,m);
  YN = fit(M,x,m);
  l=(n-1)*x/m; 
  D = YM(1:l)-YN(1:l);
  E = sqrt(sum(D.*conj(D))/l);         
endfunction
n=1024;
Z=arrayfun(@(x) error(x,x*2,n*10),[2:n]);
Y=arrayfun(@(x) 0.25/sqrt(x),[2:n]);
loglog([2:n],Z,'-',[2:n],Y,'-');

The resulting plot shows that the error in the least-squares approximation error does decrease, but nowhere near as quickly as 1/n.

EDIT: It seems to be a good fit sqrt(n)/4. I've adjusted the fit line in the program.

"Unique" band limited function (Digital Show and Tell)

Reply #9 –
I haven't looked at your code, but you might want to try something really simple (matlab code):

Code: [Select]
Fs = 192000;
M = 128;
t = (0:1:20-1)./Fs;
x = rand(1, length(t)).*2-1;
t2 = (0:1:20*M-1)./(Fs*M);
x2 = zeros(1, length(t2));
x2(1:M:end) = x;

f = sinc((-200:1:200).*0.9./M).*kaiser(401, 12)';
y = fftfilt(f, x2);

plot(t,x,'o',t2-200/(Fs*M),y)

That's simple oversampling by a factor M. The worse the filter f, the bigger the "reconstruction" error.
"I hear it when I see it."

"Unique" band limited function (Digital Show and Tell)

Reply #10 –
julf and Saratoga, thank you both so much for sharing your stories. Good to know there are people a lot more educated than myself out there but I sure wish I could do mathematics like that. Regards.
 

"Unique" band limited function (Digital Show and Tell)

Reply #11 –
Thanks for the reply, xnor.

I need to read up on these signal-processing functions you are using. I don't have a background in signal processing, so it may take me a day or two to learn the math behind FIRs. I'll reformulate my post once I understand things from your perspective.

"Unique" band limited function (Digital Show and Tell)

Reply #12 –
Quote from: terpstra on
I understand, of course, that using a DFT as I've done it is periodic. Thus my input functions already have the infinite support you desire. However, in the real world of a DAC that produces the resulting analogue function, there must be some sort of windowing going on. The DAC is not an oracle having perfect knowledge of the past and future. Thus, it must do something like I did: select a limited number of samples and fill in the gaps from there.

Sure, DAC typically uses a window of 20–60 samples.

Quote from: terpstra on
What is your 't'? Is it the number of samples?

Yes, it is time in no particular units. 1/t is the speed of decay of a sin(x)/x interpolating function that is used by both Nyquist theorem and DACs.

"Unique" band limited function (Digital Show and Tell)

Reply #13 –
Quote from: terpstra on
The resulting plot shows that the error in the least-squares approximation error does decrease, but nowhere near as quickly as 1/n.

That's because FFT does not give you a proper reconstructed function. Try using the interpolant suggested by xnor.

"Unique" band limited function (Digital Show and Tell)

Reply #14 –
Well, you can interpolate with the FFT in a quick and dirty way, but the result will be equally dirty.

Just zero pad in the frequency domain, for example double the bins to oversample by 2x, then ifft. You will get aliasing in the time domain since the FFT simply wraps around.
It's basically fast convolution with a really steep filter that has magnitude 1 up to Fs/2, and 0 above that.
"I hear it when I see it."

"Unique" band limited function (Digital Show and Tell)

Reply #15 –
So, I've finished educating myself on this topic.

The source of my confusion was that I didn't understand that the DTFT took a discrete time signal to a continuous frequency representation. Now I understand that there are 4 different ways doing a Fourier transform:
1. Continuous unbounded time <=> Continuous unbounded frequency (Fourier transform)
2. Continuous bounded time <=> Discrete unbounded frequency (Fourier series)
3. Discrete unbounded time <=> Continuous bounded frequency (Discrete time fourier transform / DTFT)
4. Discrete bounded time <=> Discrete bounded frequency (Discrete fourier transform / DFT)

I only knew about #1, 2, and 4 and I didn't really understand how they were related. Now I understand that you can transform between bounded and unbounded via periodic summation. Once both are bounded, you immediately get the DFT. The sampling theorem is just an obvious consequence of the DTFT, but I was trying to shoe-horn it into my understanding of the DFT.

The sinc(x) function just pops out of a DTFT followed by its continuous-time inverse:
x(t) = 1/2pi Int_w [Sum_m x[m] e^{-iwm)] e^{iwt} dw = Sum_m x[m] sinc(t-m)

I knew that you could fit a finite series of data points by a finite series of sine waves (DFT), so that was what I imagined the frequency domain of the DTFT to look like: a series of discrete impulses. Now that I understand it is continuous, I realize that DTFT is actually the sum of one helix spanning the frequency band per data point. Obviously, a helix turns into a sinc function when you invert the transform.

To answer my question about convergence:
Code: [Select]
function Y = fit_sinc(y,m)
  n=size(y,2);
  r=n/m;
  s=m/n;
  Z=sinc(-n+r:r:n-r);
  X=zeros(1,m);
  X(1:s:m-1)=y;     
  Y = conv(X,Z,"same");
endfunction   
n=64; 
m=4096;
X=rand(1,n);
plot([0:n-1],X,'o',[0:n/m:n-1/m],fit_sinc(X,m),'-',[0:n/m:n/2-1/m],fit_sinc(X(1:n/2),m/2),'-')


The sinc function fits nice and local due to the the division by x. My original concern that the fit is not unique still stands, despite claims to the contrary made in this thread. As long as the time series is discrete AND bounded, the DTFT does not actually provide a unique fit. The "wiggles" between samples from potential sinc functions corresponding to data points beyond the cutoff are missing. However, the error due to this truncation is bounded above by an alternating harmonic series, and quite ok in practice (due to random cancellation). Nevertheless, one could construct a worst-case signal that has an arbitrarily large reconstruction error due to truncation.

Anyway, I get it now. Yay.

EDIT: Thanks for the pointer to sinc(x), xnor. It was the realization that the Fourier Transform of sinc is continuous that started me on the golden path.

"Unique" band limited function (Digital Show and Tell)

Reply #16 –
I'm glad you get it, but...

The point of sampling theory is that you need to keep everything below fs/2, and completely remove everything above fs/2. In the frequency domain, that's a brick wall low pass filter with a cut off frequency of fs/2. To avoid changing the signal, it needs to be linear phase. That's the frequency domain representation of the filter completely defined (except at exactly fs/2). You seem to be getting a bit stuck in the maths of proving that the time domain representation of that filter is a sinc pulse, but it is. If you can't do the maths, use Fraunhofer diffraction as an analogy for a Fourier Transform (rather than the other way around) and look at the result on a piece of card from a monochromatic light source shining through a suitably small slit*. There it is  One sinc function, before your very eyes, not an equation in sight.

* your suitably small slit is equivalent to the frequency response of your low pass filter, passing everything form -fs/2 through DC up to +fs/2, and blocking everything outside of that range.

(My old maths teacher would kill me at this point.)


You're worried about truncation. The signal and the sinc function can both be as long as you need them to be. The sinc function dies away pretty quickly, so a few thousand taps gets you beyond audibility. The signal is whatever it is - don't imagine it's undefined beyond the ends - imagine it's zero. You now have a unique theoretical solution, and a unique practical one (to whatever limit you chose for the length of the sinc-based filter).

Cheers,
David.

"Unique" band limited function (Digital Show and Tell)

Reply #17 –
I'm not really sure what you're driving at. I proved in my previous post that a sinc is correct, so I certainly agree with you..? The point of the sampling theory IMO is that it states the conditions under which perfect reconstruction is possible. So it is pretty obvious to me that the right ? thing to do for digital audio is use an analogue low-pass filter to chop off the high range, with a transition centered around 36kHz or something. Then sample at 96kHz to digital (assuming your analogue IIR filter had enough attenuation above 48kHz), perform a FIR by convolution with a sinc function to precisely chop off all frequencies above 23.9kHz and then just drop every second sample to get a perfect 48kHz, with no aliasing....

Quote from: 2Bdecided on
You're worried about truncation. The sinc function dies away pretty quickly, so a few thousand taps gets you beyond audibility.

It only dies away quickly because the error from the sinc functions cancel out randomly. Imagine your signal beyond the point of truncation was 1, 0, 1, 0, 1, 0, ... Then the error between your last two samples is: ~ 1/2 + 1/4 + 1/6 + 1/8 + 1/10 + ... That is a harmonic series and it diverges. That's all I was moaning about.

Quote
The signal is whatever it is - don't imagine it's undefined beyond the ends - imagine it's zero.

Well, it is not zero. It is undefined. I agree that in practice truncation won't matter because the errors are not going to be aligned.

"Unique" band limited function (Digital Show and Tell)

Reply #18 –
Quote from: terpstra on
It only dies away quickly because the error from the sinc functions cancel out randomly. Imagine your signal beyond the point of truncation was 1, 0, 1, 0, 1, 0, ... Then the error between your last two samples is: ~ 1/2 + 1/4 + 1/6 + 1/8 + 1/10 + ... That is a harmonic series and it diverges. That's all I was moaning about.
In a way, that particular behaviour is also as expected. Your signal is at (what will be) Nyquist. The response exactly at Nyquist is undefined, and signals at Nyquist are not recoverable.

Quote from: 2Bdecided link=msg=0 date=
Quote
The signal is whatever it is - don't imagine it's undefined beyond the ends - imagine it's zero.

Well, it is not zero. It is undefined. I agree that in practice truncation won't matter because the errors are not going to be aligned.
Let us accept that at least some transports and DACs will treat the digital silence before you hit play on the CD player as being a string of 0s. There you go. It's zero

I know this isn't helping your equations. Someone else can do that.

Cheers,
David.