One other topic to keep in mind is how you are performing the DFT (or FFT) and subsequently how you process the signal. Different techniques result in different assumptions about the signal, which can significantly effect the result, especially near the Nyquist frequency. For instance, a few years ago, I was looking at the Hilbert transform of a signal to obtain it's envelope. Using the standard MATLAB frequency sampling based Hilbert transform resulted in significant artifacts, especially at the beginning and end of the signals. I instead used a Hilbert filter based on a Remez Exchange algorithm, which designs the filter by approximating the desired frequency response at all frequencies below the Nyquist criteron (not just the sampled ones). This significantly reduced the artifacts observed. You may find that a similar concept will be useful in this application.

You can try sinc interpolation. When used with enough samples values it should work fine.

If data are only around f=1/2 you can shift signal spectrum, so the data are located at lower frequency (for example multiplying signal by (-1)^n) and then use any interpolation technique you like with better results. But getting signal back at correct frequency if it is real-valued might be tricky.

In general case it might be better to carefully design interpolation filter since in your case imple sinc interpolator might require extremely long impulse response. First you need to decide on inrepolation filter specification - passband and stopband location. It is important to carefully select transition band since its width determines impulse response length. With interpolation filter designed you can simply use zeroinserting followed by filtration (or polyphase implementation of such approach). If expected impulse response length is to long for the design procedures, you can use Interpolated-FIR or Frequency Masking techniques.

If the resampling ratio is irrational then I would recommend fractional delay filters, though the filters should be carefully designed to obtain overall filter fullfilling your requirements.

This is a rather open question. It all depends on the signal you are sampling. If you look at the problem from a simple Nyquist view, i.e. the signal is completely contained within the bandwidth (f), and the sampling signal is a pure impulse, then the signal should be recoverable 'perfectly'. However, in reality both these are not true, in which case you need to take into account the target signal frequency response, the sampling signal frequency response, and only then you will be able to determine the best way forward.

As an example, consider the case where the target signal is mostly contained within the bandwidth (f) although beyond the frequency (f) the signal power is low but not 'zero'. This non-zero signal will be aliased back into the 1st Nyquist zone which will contribute to interference. Also combined with this with also assuming that your sampling signal is not an impulse, but non-perfect square wave with a 'mostly' sinx/x profile but not exactly. The resulting sampled signal will the product of the two signal frequency responses, (i.e. aliased target signal * sampling signal frequency response).

You can imagine that he result will not be a perfectly Nyquist sampled result. If on top of this you start adding implementation issues, such as digital sampling, quantisation issues, noise, non-linearity's etc then you move even further away from the ideal.

With regards to your choice of the cubic spline, it is a reasonable choice, but keep in mind that this filter also has a non-flat frequency response which you are invariable convolving with your input signal. If you analyse this cubic spline response you may recognise why your output sampled signal is not what you expected. In essence it will be applying low pass filtering of your target signal with a filter profile which is not a perfect 'brick wall'.

I would do an impulse analysis of your interpolator and try to understand my target signal a little better in order to determine why the output is what it is. You may very well notice by looking at the spectrum of your output signal that the interpolator is contributing to a significant distortion which will explain your observed output.

Hope this helps.

Christopher.

You cannot, of course, "recover the original analog data".

I suppose that the question is: how do you obtain an equivalent series of digital data at a different (greater) sampling rate.

I would do FFT of tapered data interval, then zero padding in Fourier space, inverse FFT, recovering for tapering (i.e. multiplication by reciprocal of the tapering function).

Depending on how much your signal is near to the Nyquist rate, you will need to scratch some amount of data at both ends of your time interval, and compensate by overlapping next interval (overlap-save method).

I obviously assume that there was no aliasing problem in the first place when the signal was sampled (good anti aliasing very near to the Nyquist frequency imply very narrow skirt).

This may be a silly question by me, but can I ask about the spectrum?

When you say "close to the Nyquist rate" can I confirm that you mean that the spectrum is relatively full, -fN ..fN, & flat, white or is the data you want to restore just near Nyquist?

I guess as Marek Blok has suggested, resampling, simple spectral inversion, or similar tricks may be useful.

Otherwise, in general, the longer the interpolant, the closer to Nyquist you can aspire to get. I guess your limit may have been determined by the original anti-alias filter, or equivalent process) prior to sampling?

I agree with L. Baggio. The simplest solution is to go to frequency domain, add as many zeros as you need and finaly go back to time domain. This way you apply ideal (with square shape) lowpass filtering and also you increase sampling freq. So, in infinity you can obtain analog signal ;)

The other problem is with the digital data 'sampled at Nyquist freq' as you mentioned. If the original analog signal before sampling have some - even small - portion of the energy for freq bigger or equal to sampling freq you obtain aliased digital data and to the best of my knowledge there are no method to resolve this problem.

Regards,

Boguslaw

Hi. Thank you all for the tips. My digital signal is wideband and I am intererested in extracting the energetic amplitude of all the spectral bands. The only problematic components are those contained in the highest bands, next to the Nyquist rate. Lt's suppose we have a pure tone of frequency almost 1/2, then its envelope will not be constant as expected. Therefore, I propose to obtain a "more analog" version of this band. I think zero padding in the spectrum will work, since oversampling is not enough.

% simple example for Matlab

N=1024; % original data length

x=normrnd(0,1,1,N); % original signal with the spectrum from 0 to Fs / 2

% where Fs sampling

figure(1); periodogram(x) % spectrum of the original signal

% Time domain

y=zeros(1,N*4); % upsampling (add zeros)

for i=1:N y(i*4-3)=x(i); end

figure(2); periodogram(y)

B=fir1(64,0.25); % design interpolaton filter

freqz(B,1) % Digital filter frequency response

z=filter(B,1,y); % filters the data

figure(3); periodogram(z) % spectrum interpolated signal

% frequency domain

xf=fft(x);

figure(4); plot(abs(xf));

yf=zeros(1,4*N); % add zeros in frequency domain

for i=1:512 yf(i)=xf(i); yf(4097-i)=xf(1025-i); end

figure(5); plot(abs(yf));

zt=real(ifft(yf));

figure(6); plot(1:4096,z,1:4096,zt); % compare interpolation signals

figure(7); periodogram(zt);

I would use the "sinc" interpolation with most or all of the signal being considered (i.e. inside the window). The sinc interpolation is more or less the definition of the signal theorem and is easy to apply (in my opinion).

Edit: It should be mentioned that since sinc is equivalent to the definition of the discrete sampling theorem, every other function will add frequencies in comparison to sinc. Whether you want those or not you have to decide.

Take a step back to see what you are actually doing: You have a signal that is sampled and you try to recover the original. To do that you define a base set of functions and then you choose for every local sample the best fitting parameters in your base set. This set of base functions should fit well to the original signals to get the best results.

From your test with the splines functions set it seems that they do not really fit the form of the original signals. So the thing to do would be to write down what you know of the original signals form, find a functional base that fit them well and then try to fit the sampled data with this new base.

You may want to have a look at literature in the field of 'compressive sensing'. Your problem is one of the main topics there.

Hello You can use Spectral Interpolation which is basically adding enough zero to "compress" the signal spectrum away from the Nyquist rate. I have used this method for irrational sample rate adjustment and works perfect. One should notice that we are not adding any information; the compression happens because with the added zeros(at one end of the signal spectrum), the interval -pi to pi is assigned to a larger number of points.

Masoud

Adding zeros in DFT domain is fast approach and does not change the information but this doesn't mean that this is a good interpolation technique. This method has two drawbacks: (1) it results in time domain aliasing and (2) actual frequency response of the inrepolation filter suffers from Gibs phenomenon. The (1) can be reduced a bit by using large zeropadding at the begining in the time domain. Time domain aliasing will be smaller but you need remember that the impulse response the interpolation filter is the aliased sinc of the same length as the DFT length. Therefore with circular convolution we obtain large transitional states. The (2) can be improved with careful selection of addtional samples of the frequency response of the interpolation filters in the transition band just like in classic frequency sampling filter design technique. This would also help with (1) since the resulting impulse response would decay faster at the ends resulting in shorter efficient impulse response length. Generally this technique is simple and fast but it is not precise. If a designer needs more control other the performance of the interpolator, I would rather advice the use of optimal interpolation filters for example designed using Parks-McClellan method and use classic zeroinserting - filtering - decimation approach. If one prefers frequency domain implemantation I sugest "Convolution-Based Trigonometric Interpolation of Band-Limited Signals" by J. Selva. Although it is described for bandlimited signals, it can be actually used for fulband signals with distorsions located in transition band. For nearly nonoversampled signal you just need to use filter with passband allmost up to f = 1/2 (depending on your needs).

What can be achieved in general in terms of reconstruction of a non-ideally sampled signal has been thoroughly investigated by Michael Unser at EPFL and his collaborators.

See http://bigwww.epfl.ch/publications/;

for starters I recommend to have a look at

M. Unser, "Sampling—50 Years After Shannon," Proceedings of the IEEE, vol. 88, no. 4, pp. 569-587, April 2000.

and then have a look at the very recent papers on sampling and interpolation from the same group.

You might also have a look at my group's work on optimal interpolation, reconstruction from non-ideal sampling:

http://www.vsi.cs.uni-frankfurt.de/research/current-projects/optimal-interpolation-reconstruction-from-non-ideal-sampling/

The Shannon conditions (precisely bandlimited processes and ideal Dirac-sampling) are highly relevant theoretical special cases, but _real_ signals are not truly bandlimited, and most processes will decay in power density already significantly way before the Nyquist frequency. It is reasonable to use this knowledge in order to design a reconstruction scheme which does not give you a perfect reconstruction (this is only possible in Shannon's perfeect world), but the optimum one, e.g. under a least squares criterion.

Agree with L Baggio, but thanks for Mester update

One other topic to keep in mind is how you are performing the DFT (or FFT) and subsequently how you process the signal. Different techniques result in different assumptions about the signal, which can significantly effect the result, especially near the Nyquist frequency. For instance, a few years ago, I was looking at the Hilbert transform of a signal to obtain it's envelope. Using the standard MATLAB frequency sampling based Hilbert transform resulted in significant artifacts, especially at the beginning and end of the signals. I instead used a Hilbert filter based on a Remez Exchange algorithm, which designs the filter by approximating the desired frequency response at all frequencies below the Nyquist criteron (not just the sampled ones). This significantly reduced the artifacts observed. You may find that a similar concept will be useful in this application.

If you look at the autocorrelation function, you may find that you have what you need