Thursday, June 26, 2008

sinc interpolation in MATLAB

UPDATE: The examples given here are meant to give mathematical insight into how sinc interpolation works by using a finite-time APPROXIMATION. Sinc interpolation necessarily involves a circular convolution, which is not a finite computation in the time domain. If you actually need to do sinc interpolation, use the interpft function. It does an FFT, pads the FFT with zeros, and does an IFFT. Consequently, it is VERY FAST. Moreover, using an FFT (or DFT, in general) is the only way to use a finite computation to do a sinc interpolation (or circular convolution in general). Just make sure that the resampled time vector you use has a length that is an integer multiple of your original time vector. Also make sure that it lands on the same points as your original time vector (i.e., it should only add new points between old points).

I've also created a interpftw that does the same job as interpft but allows you to reconstruct samples from different aliasing windows. Whereas interpft pads the end of the FFT with zeros, interpftw pads both the beginning and the end. In other words, the former is a low-pass filter and the latter is a bandpass filter.
If you search Google for sinc interpolation in MATLAB, many pages will reference the sinc_interp example from John Loomis. Unfortunately, I've found few sites that recognize that the function is not meant to do general purpose sinc interpolation. That is, it makes a few assumptions about the sampling rates that may not be evident to the average user.

So, I'm giving some of my students this example:
% Ideally "resamples" x vector from s to u by sinc interpolation
function y = sinc_interp(x,s,u)
    % Interpolates x sampled sampled at "s" instants
    % Output y is sampled at "u" instants ("u" for "upsampled")


    % Find the sampling period of the undersampled signal
    T = s(2)-s(1);

    for i=1:length(u)
        y( i ) = sum( x .* sinc( (1/T)*(u(i) - s) ) );
    end

    % Make sure y is same shape as u (row->row, col->col)
    y = reshape(y, size(u));
end
Here's a vectorized (i.e., MUCH FASTER) version:
% Ideally "resamples" x vector from s to u by sinc interpolation
function y = sinc_interp(x,s,u)
    % Interpolates x sampled at "s" uniformly spaced instants
    % Output y is sampled at "u" uniformly spaced instants
    % ("s" for "sampled" and "u" for "upsampled")
    % (consequently, length(x)=length(s))


    % Find the period of the undersampled signal
    T = s(2)-s(1);

    % The entries of this matrix are each u-s permutation.
    % It will be used to generate the sinc transform that will
    % be convolved below with the input signal to do the
    % interpolation.
    %
    % (recall that u(:) will be a column vector regardless
    % of the row-ness of u. So u(:) is a row, and s(:) is a
    % column)

    sincM = repmat( u(:), 1, length(s) ) ...
           - repmat( s(:)', length(u), 1 );

    % * Sinc is the inverse Fourier transform of the boxcar in
    % the frequency domain that was used to filter out the
    % ambiguous copies of the signal generated from sampling.
    % * That sinc, which is now sampled at length(u) instants,
    % is convolved with the input signal becuse the boxcar was
    % multipled with its Fourier transform.
    % So this multiplication (which is a matrix transformation
    % of the input vector x) is an implementation of a
    % convolution.
    % (reshape is used to ensure y has same shape as upsampled u)

    y = reshape( sinc( sincM/T )*x(:) , size(u) );
end
My function sinc_interp resamples the data in the x vector. The original time vector is given by the s vector and the new time vector is given by the u vector.

Alternatively, you can try this more advanced one that lets you pick which aliasing window you want to use to reconstruct (i.e., it's an ideal bandpass filter rather than just a low-pass filter).
% Ideally "resamples" x vector from s to u by sinc interpolation
function y = sinc_interp(x,s,u,N)
    % Interpolates x sampled sampled at "s" instants
    % Output y is sampled at "u" instants ("u" for "upsampled")
    % Optionally, uses the Nth sampling window where N=0 is DC
    %     (so non-baseband signals have N = 1,2,3,...)


    if nargin < 4
        N = 0;
    end

    % Find the period of the undersampled signal
    T = s(2)-s(1);

    for i=1:length(u)
        y( i ) = ...
            sum( x .* ...
                ( (N+1)*sinc( ((N+1)/T)*(u(i) - s) ) - ...
                  N*sinc( (N/T)*(u(i) - s) ) ) );
    end
end
Here's a vectorized (i.e., MUCH FASTER) version:
% Ideally "resamples" x vector from s to u by sinc interpolation
function y = sinc_interp(x,s,u,N)
    % Interpolates x sampled sampled at "s" instants
    % Output y is sampled at "u" instants ("u" for "upsampled")
    % Optionally, uses the Nth sampling window where N=0 is DC
    % (so non-baseband signals have N = 1,2,3,...)


    if nargin < 4
        N = 0;
    end

    % Find the period of the undersampled signal
    T = s(2)-s(1);

    % When generating this matrix, remember that "s" and "u" are
    % passed as ROW vectors and "y" is expected to also be a ROW
    % vector. If everything were column vectors, we'd do.
    %
    % sincM = repmat( u, 1, length(s) ) - repmat( s', length(u), 1 );
    %
    % So that the matrix would be longer than it is wide.
    % Here, we generate the transpose of that matrix.

    sincM = repmat( u, length(s), 1 ) - repmat( s', 1, length(u) );

    % Equivalent to column vector math:
    % y = sinc( sincM'(N+1)/T )*x';

    y = x*( (N+1)*sinc( sincM*(N+1)/T ) - N*sinc( sincM*N/T ) );
end

10 comments:

Royi said...

Do you have a proof why is Zero Padding in the Frequency domain is Sinc interpolation.
Intuitively it's easy to understand, yet I'm looking for a Mathematical proof.

If you let me know by email it would be great:
royia@walla.co.il

Thanks.

Ted said...

Royi – such a "proof" exists in most digital signal processing (DSP) textbooks and probably many textbooks on (discrete-time) systems theory in general. In the end, the relationship is "because" multiplication by a boxcar in the frequency domain is equivalent to convolution with a sinc function in the time domain. Because interpolation reduces the sampling frequency, it introduces ambiguity about frequency content because the signal is now "allowed" to change more quickly between the old points. Sinc interpolation picks those intermediate points to be consistent with boxcars (i.e., idealized low-pass or bandpass filters) in the frequency domain.

If you start at the frequency-domain series and write out the expression for the inverse Fourier transformation after padding the series with zeros at both ends, you will get the expression for sinc interpolation directly.

Perhaps you are looking to "prove" something else. For example, you might want to prove that the resulting time-domain signal matches the original time series at all of the original times and places data in between along smooth curves that connect adjacent datum. These sort of proofs/demonstrations can also be found in most DSP textbooks.

Some things that may ease your mind...

1. Decreasing sampling time necessarily widens the Nyquist frequency. That is, decreasing sampling time increases the possible unambiguous bandwidth of signals being represented.

2. Interpolation should not ADD any frequency content. That is, no information should be added by interpolation. Instead, redundant information should be "filled in" between points.

3. To maintain the original frequency domain, I should put an ideal low-pass filter around it that has a brick wall edge at the old Nyquist frequency. That way I am sure I do not add any energy at the higher frequencies that I have introduced by decreasing my sampling frequency.

4. It is this idealized (i.e., brick wall) implementation of a low-pass filter that results in zero padding. Everything in the original frequency-domain window remains the same, and anything in the higher frequencies is "filtered out" and replaced with zero.

You can extend this explanation to bandpass filters as well. That is, if you assume your original frequency-domain content is not at the baseband (i.e., not centered at DC), then you want to make sure your wider frequency window doesn't include any baseband nor high frequency energy. Hence, you put your extra zeros not only at high frequency but also near DC (baseband).

I hope that sets you on the right track to getting your questions answered.

Royi said...

Hello.
I couldn't find it in any DSP book.
Most of them doesn't even say how to zero pad correctly in the frequency domain.

I want to see a rigorous proof how zero pad is the equivalent of higher sampling frequency (Assuming Limited band function).
Just as you described it.

Could you point me for such a proof?

Thanks.

Ted said...

Royi - I am really surprised that you couldn't find the information that you need in DSP-specific textbooks you have. This is fairly simple undergraduate level material typically found in an introductory course on the subject. Unfortunately, most of my signal processing texts are locked up in my office, and I rarely go there nowadays, and so I don't have a particular reference to give you.

I think I said this before, but perhaps the problem is how you are framing your question. When you ask for a "proof", you have to come up with a statement that you want to prove (e.g., "if f is differentiable, then f is continuous"). I'm not quite sure I see that kind of statement here (at least not with a requisite amount of detail). Again, I stress that you should WRITE OUT the DEFINITION of the inverse (discrete-time) Fourier transform and ADD ZEROS onto the end of your frequency domain data. It will be simplest if you use the complex exponential form, but you will need to be comfortable handling complex numbers. Otherwise, you will have to be comfortable with trigonometric identities. Either way, with a little algebra, you should see that you get the same sampled signal but at a higher sampling frequency. It comes right out of the definition of the inverse transform (just remember to repeat your frequency domain window appropriately when you're doing your sums). I have a feeling that you will only be satisfied by working out the math on paper yourself. You will likely not find this exact "proof" in a book because it is basically an exercise in algebra.

I did a Google Book search and got some promising results. Try searching on books.google.com for "zero padding frequency domain". If I recall correctly, the book

_A course in digital signal processing_ by Boaz Porat

was my undergraduate testbook in DSP. As you can see, on page 106 (section 4.5), there is a section on zero padding in the frequency domain that I think might help you. Unfortunately, Google books only offers a snippet view. However, you should be able to find this book in any university library.

You should also get a result from the Google Book search:

_Digital Signal Processing using MATLAB and Wavelets_ by Michael Weeks

Look on page 202 in section 6.4 ("Zero padding"). This section discusses zero padding in the time domain to provide a better frequency resolution. However, there is a duality between the two domains. Hence, you can switch "frequency" and "time" and find an argument for how padding in the frequency domain adds no information but increases time-domain resolution. The author extends the argument to discrete time without any mathematical details, but you can use his continuous time example to build the case in discrete time.

Another Google Books result:

_Signal processing: principles and implementation_ by S. V. Narasimhan, S. Veena

gives discrete-time details on page 75 (section 4.9: Zero padding).

I'm not sure these references will have exactly what you need, but hopefully they will help you when you are constructing your own proof (i.e., writing out the IDFT of zero-padding signals and convincing yourself that you get a sampled time-domain signal matching a (sinc-)interpolated version of the original). Good luck!

Royi said...

Well, let's define what I'm asking.

Padding in Frequency Domain is equivalent of higher frequency of sampling, meaning Interpolation.
This is easy to prove.

What's hard to prove (Rigorously, not by intuition) is how should this padding be done.
Let's talk about the even case.
We should take the center item (Of the DFT) split it into half and zero pad in between.

Could you locate a proof that goes all the way by showing that Padding is Interpolation and conclude you need to do it by dividing the center item into 2?


So what I'm asking is a proof to show the right method for doing the padding in the frequency domain.

I once did it, not so rigorously.
The main idea is using Sinc as interpolator.
In the DTFT domain higher frqeuncy will result by a narrower bandwidth (Relative to the normalized frequency Pi). Multiplying it by "Window". Yet one should remember that "Window" Fourier Series converge to half the amplitude at the edges, hence the dividing.

It's not Rigorous hence I'm looking for something deep.

I hope this time I'm clear.

Thanks.

Ted said...

Royi --

1. You cannot take for granted that "padding" in the frequency domain is equivalent to interpolation in the time domain. It is true that because the transformation is an operator, if you have N samples in one domain, you will have N samples in the other (by definition). But that's all you can say about that. Depending on where you put your zeros, you can greatly change the way the time-domain signal looks.

2. "Interpolation" is a general term. There are lots of different ways to place points between other points. "Sinc interpolation" is one method. There is also "linear interpolation" (of the "first-order hold" variety), "zero-order hold", etc. Zero-padding just provides a way to find points between other points so that no additional harmonics are added by the interpolation process.


Here is what I recommend you do... ("Step 2" below is probably closest to what you want)

Step 1 - From Time to Frequency:

(a) Define a generic continuous-time signal f(t) that is band-limited to a window from -1/(2T) to 1/(2T) (e.g., where "T" is in seconds)

(b) Sample f(t) every T seconds, giving you the discrete-time signal f_1[k] := f(kT)

(c) Additionally, sample f(t) every T/n seconds, giving you the discrete-time signal f_n[k] := f(kT/n)

(d) Use the definition of the DFT (i.e., the finite sum of complex exponentials) to write out the DFT of (b) and (c).

Because f(t) is band limited, the DFT of (c) will have zeros at the discrete frequency points corresponding to frequencies above 1/(2T). Hence, the "right" interpolation from f_1 to f_n leads to zero-padding to the right of the DFT.


Step 2 - From frequency to time:

(a) Start with a [-1/(2T),1/(2T)) band-limited DFT F[k] and takes its IDFT to generate a time-domain f[k]

(b) Generate the DFT F_n[k] by zero-padding the right-hand side of F[k] out through [-n/(2T),n/(2T)) and take the IDFT to generate a time-domain f_n[k]

(c) Verify that f[m]=f_n[m*n]

By definition, f_n[] is an "interpolation" of f[], and both signals remain band-limited between -1/(2T) and 1/(2T). Additionally, you can verify that each f_n[k] is a sum of the original f[k] values with weightings that can be written as sinc functions. Hence, we have "sinc interpolation."

Contrast "sinc interpolation" with (*) the "linear interpolation" of a "first-order hold" circuit, or (**) a "zero-order hold" circuit. In both cases, these hold mechanisms introduce harmonic distortion. In time-domain terms, these hold mechanisms introduce sharp edges. You don't get those "sharp edges" with sinc interpolation because sinc interpolation is naturally a low-pass filter.

You can extend these exercises to the bandpass case as well.

Anonymous said...

"If you actually need to do sinc interpolation, use the interpft function"

But those aren't the same thing. interpft() assumes that the signal is periodic and repeats outside the known values. Your sinc interpolation assumes that everything outside of the known values is 0. If you pad the signal with zeros it will approximate the true sinc interpolation more and more, but I don't think it is ever perfect. Is there a way to make a perfect reconstruction with FFTs?

Also your link is broken. http://www.mathworks.com/help/techdoc/ref/interpft.html

Ted said...

Anonymous -- Thanks for catching the dead links. I've updated them.

Regarding your other comments, keep in mind that the padding is done in the FREQUENCY domain and not the time domain. With more padding at the edges of the frequency domain, intermediate points are added in the time domain that are consistent with the given spectrum.

Sinc interpolation naturally assumes that the signal is periodic outside of the given interval. In particular, sinc interpolation assumes that the signal is bandlimited and "fills in" the blanks consisted with the information captured in the known frequency band. If it is bandlimited, it has finite support in the frequency domain and thus cannot have finite support in the time domain.

On the other hand, if you assume a time-limited signal (i.e., a signal that has finite support in the time domain), you do have to pad the signal with an infinite number of zeros in the time domain. The result interpolates in the frequency domain. Frequency components are added in between *and around* the original components that are still consistent with the ever-padded time-domain signal. The limiting curve is the smooth Fourier transform for the time-limited signal. That is, padding with zeros in the TIME domain approximates the Fourier transform in the frequency domain.

Jonathan Dandy said...

Thanks for the post - I found it useful. I think you have an error in the comments of the vectorized low-pass code.
"% (recall that u(:) will be a column vector regardless
% of the row-ness of u. So u(:) is a row, and s(:) is a
% column)"
should be
"% (recall that u(:) will be a column vector regardless
% of the row-ness of u. So u(:) is a column, and s(:)' is a
% row)"

savs95 said...

You said to make sure that when using interpft function, the interpolated points lands on the sampled points given before. How to do that, as I can't get that through my code.

My code is here
f1 = 100;
fs = 200;

time = linspace(0, 3/f1, 200);
x = cos(2*pi*f1*time);
t = linspace(0,3/f1,10);
y = cos(2*pi*f1*t);

vq1 = interpft(y,200);
plot(t,y,'ro',time, vq1);


Here it's not landing on same point