[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: data resampling
Hi,
to understand the difference between the three methods, looking at a
swept sine could also be helpful. The main difference between decimate
and resample in this case is the type of filtering applied before the
downsampling. These filters can be adjusted. For example, an order 4
decimate can be more similar to the results of resample with default
settings than to decimate with default settings. For my applications,
resample with N set to 50-100 is usually fine.
fs = 8000;
len_sec = 3;
t = [0:1/fs:len_sec];
y = sin(2*pi*( fs/(4*len_sec) * t.^2 ));
lenwin = 2^nextpow2(fs * 0.04);
window = hann(lenwin);
[z,f,t] = spectrogram(y,window,round(lenwin*3/4),[],fs);
figure; imagesc(t,f,abs(z)); title('input')
fs = fs/2;
lenwin = 2^nextpow2(fs * 0.04);
window = hann(lenwin);
[z,f,t] = spectrogram(decimate(y,2),window,round(lenwin*3/4),[],fs);
figure; imagesc(t,f,abs(z)); title('decimate default')
[z,f,t] = spectrogram(decimate(y,2,5),window,round(lenwin*3/4),[],fs);
figure; imagesc(t,f,abs(z)); title('decimate order 5')
[z,f,t] = spectrogram(resample(y,1,2),window,round(lenwin*3/4),[],fs);
figure; imagesc(t,f,abs(z)); title('resample default')
[z,f,t] = spectrogram(resample(y,1,2,50),window,round(lenwin*3/4),[],fs);
figure; imagesc(t,f,abs(z)); title('resample length 50')
[z,f,t] = spectrogram(downsample(y,2),window,round(lenwin*3/4),[],fs);
figure; imagesc(t,f,abs(z)); title('downsample')
Sebastian
On 15/08/2014 08:26, Alain de Cheveigne wrote:
Hi Hari,
Note that matlab's resample() has inadequate antialiasing. Try:
figure(1); clf;
x=randn(10000,1);
pwelch(decimate(x,2)); hold on
pwelch(resample(x,1,2));
c=get(gca,'children'); set(c(1), 'color', [0 0 0]);
legend('decimate','resample');
decimate() does the right thing which is to apply a lowpass filter with a sufficiently low cutoff so that power is negligible at Nyquist frequency and beyond. resample()'s apparently applies a filter that attenuates by 3dB at Nyquist, so that there may be serious aliasing if the signal is wideband.
I filed a bug report to MathWorks, they answered that they think it's OK (!). Annoyingly decimate() works only for single channel data.
It's an example of matlab's tendency to have lots of different functions for a similar purpose, each limited or buggy in a different way (see http://abandonmatlab.wordpress.com/2012/08/07/matlab-cant-read-plain-text-data-out-of-a-wet-paper-bag/ for an enlighening rant).
Concerning the original post, the best form of downsampling depends on your requirements. If the goal is to preserve time-domain properties of the waveform, it may be best to simply replace sets of N samples by their average (where N is the downsampling factor). If the goal is to preserve frequency-domain properties (within the band below Nyquist), decimate is best. Both are similar in that each new sample is a function of a set (window) of samples of the original series. In the first case this set is a compact as makes sense (N). In the second it is much more extended, so temporal features may be blurred. In the example given, the window includes samples before and after the beginning and end of the signal (treated as zero) so the outcome is contaminated by the response to the onset and offset step. This effect is severe because of the large DC offset.
downsample() is unlikely to be useful in any situation as it ignores all but 1/N samples of the original data. It's a rather silly function because its functionality is easy to get by writing x(1:N:end).
Alain
On 14 Aug 2014, at 15:41, Hari Bharadwaj <hari@xxxxxxxxxxxxxxxxxxx> wrote:
Hi J,
downsample(.) just picks out every nth sample (300th in your case). Though visually it may look like it worked well, it is quite likely that some of the high frequency portions of the signal are now aliased into the low frequency portions.
resample(.) and decimate(.) guard against this aliasing by filtering first.. But filtering would cause artifacts at the edges of your signal as you don't have the input for infinite past (or in your FIR case, enough samples that you don't care about, for your filter transients to die out) . However the middle portions should be properly filtered and decimated. So if you throw out the first few samples (technically, the length of the filter used) and the last few sample, the rest should match reasonably with what you have in mind. This is not a problem with the software implementations but rather is a fundamental limitation of filtering and time-frequency uncertainty. If you care about those first few seconds, your input should include additional/extra samples at the beginning and end that you don't care about that would absorb the filter transients.
I would recommend consulting a signal processing text book for more elaborations on these issues (I.e, aliasing, filter design, resampling etc.).
Briefly, these might help:
http://en.m.wikipedia.org/wiki/Anti-aliasing_filter
https://ccrma.stanford.edu/~jos/fp/Transient_Response_Steady_State.html
The help for resample(.) talks about what kind of filter is employed: http://www.mathworks.com/help/signal/ref/resample.html
Hari
On Aug 13, 2014, at 7:05 PM, Juan Wu <wujuan22@xxxxxxxxx> wrote:
List experts,
I am attempting to down sample my data from 300 Hz to 1 Hz. I just tried two functions of Matlab: 1) a = resample(Data,1,300); 2) b = downsample(Data,300). The results are quite different between each others.
<image.png>
Apparently, the result from "downsampled" is close to the the original data. However, quite a few persons suggested using the re-sampling method - I get this information from google searching, and very agree with this view. Personally I think the downsample method is too simple and very arbitrary. I also do not believe the "Nearest Neighbor" method. I assume that the method resampling takes both "FIR interpolator and decimator" implementation - this is what I expected. Am I right? But now my output from the re-sampling method is really very terrible. So I assume that the resample function from Matlab does not do a good job for this. I am not sure whether I need change to use other softwares or try other functions in the Matlab.
Any opinions or references are much appreciated.
J
The information in this e-mail is intended only for the person to whom it is
addressed. If you believe this e-mail was sent to you in error and the e-mail
contains patient information, please contact the Partners Compliance HelpLine at
http://www.partners.org/complianceline . If the e-mail was sent to you in error
but does not contain patient information, please contact the sender and properly
dispose of the e-mail.
--
Dr. Sebastian Ewert
Centre for Digital Music, Queen Mary University of London
s.ewert@xxxxxxxxxx
Phone: +44 207 882 8287