how do I re-sample a continuous function; but maintain shape?

I am trying to resample a vector that represents a voltage time-series from 200kHz to 50 kHz using the resample function (the two vectors are shown below).

The code I've used to resample is as follows:

x = resample(x,5e4,20e4);

When I measure what we call the 'half width' (the duration of time at half the height of the voltage peak) I notice that there is an increase, essentially the half width becomes longer. The increase also varies between waveforms.

I am curious if there is anyway to reduce this discrepancy, or at least perhaps have it consistent and based on the slower sampling rate? I've tried using filter coefficients that mimick the shape of these waveforms; but that seemed to make things worse.

Thank you in advance!

Eric

[Daniel M]

Just curious what the relative percent error is? abs(width50 - width200)/width200 * 100

Note that the filtering within the resample function assumes that the signal is zero at the beginning at end points. It would help with your issue of edge effects if you subtracted out the mean before resampling and then added it back later.

[Eric Kuebler]

I plugged in your above formula and it looks like it varies between 20 and 80%, more often leaning toward 20%.

Thank you,

Eric

[Daniel M]

Wow. That's a large number. Are you sure the second version is capturing the correct half width? The figures themselves don't even look visually different. Guessing visually I would have said 1-2% difference.

How are you calculating half width?

Also note that in general when downsampling a signal, there is a loss of what is called the Q-factor of a sharp signal. (That is, a peak will get smaller in magnitude and broader in width).

[Eric Kuebler]

My apologies, I tried downsampling to 10kHz and this is where the error is much much larger.

Comparing the waveforms above (so 200kHz to 50kHz) there is a 7% difference.

[Eric Kuebler]

I am computing the half width by the following

dt = 1000/5e4;

[peakV,peakTime] = max(x);                                  % peak
[troughV,troughTime] = min(x(peakTime:end));                % trough
troughTime = peakTime+troughTime-1;                         % adjusted trough time
halfHeight = (peakV-troughV)/2;                             % half height
halfHeight = peakV-halfHeight;                              % adjusted half height
halfHeightTpre = find(x(1:peakTime)<halfHeight,1,'last');   % half height pre-peak
halfHeightTpost = peakTime+find(x(peakTime:troughTime)<halfHeight,1,'first'); % half height post-peak
hw = (halfHeightTpost-halfHeightTpre)*dt;                   % half-width

Answer 1

[Daniel M]

I see several potential reasons why this is occuring.

1) There will be general broadening of the signal due to downsampling.

2) Since the magnitude of the downsampled signal is potentially lower, then half of the max could also be lower. So you could be measuring a wider part of the signal.

3) You're measuring half-width from the left side. This combined with issue 2 could lead to measuring something like the blue line in "Capture" (attachment). I would have defined the half-width as something more like the red line, which is the HWHM from the left side.

4) You're using discrete values to calculate width. Probably not an issue because the sampling rates are still quite high, but could lead to something occuring like in "Capture2" attachment. If you find this is a big issue, then you could interpolate between the first point <halfHeight and first point >halfHeight.

I'm not really sure what to suggest other than to try using the functions pulsewidth and/or statelevels in the signal processing toolbox.

[Eric Kuebler]

1. I agree, this is why I mentioned that perhaps there is a way to do this that the difference would be based on the slower sampling rate.
2. I agree, if this is the case I would expect the same amount of error in a trace that was originally sampled at 50kHz.
3. It may be useful for us to try taking this measure from the left side (i.e., threshold).
4. This could also be the case; but agree sampling rates are still quite high.I will assess whether this is an issue.

Can you say a bit more about how I might use the function statelevels?

BTW, I didn't mean to accept this answer, I was hoping to leave it open; but hit that button by accident, good thing they put the 'Accept' button so close to the attachments!!

Answer 2

[Daniel M]

You can unaccept the answer if you want to keep it open, but also consider this answer.

(For 3, I meant you're measuring from right side, and could try from left. I think you understood what I meant).

Here are some ideas. I'm not saying they are robust, but this is certainly a potential direction to go in.

% create sample data
fs = 4e6;
t = (0:99).*1/fs;
y = [ones(1,45) 1+sin(2*pi*320000*t(1:10)) zeros(1,45)+2e4.*t(1:45) ] + 0.1.*rand(1,length(t));
% that looks ridiculous I know..., but it isn't a bad approximation to your signal

% get the statelevels
levs = statelevels(y); 
% if you call it without the output, it will display the attached figure "levels"

% get the overshoot value, need to flip y 
% I'm sure using max or findpeaks is sufficient here too. 
% I'm just showing you the utility of these functions

[os,oslev,osinst] = overshoot(fliplr(y),t,'StateLevels',levs);
% if you plot without outputs, get attached figure "overshoot"
% we care about 'upper state' and 'post-overshoot'

% get the pulse width, the crossings, and the midlevel
[width,init,fin,midl] = pulsewidth(y,t,'StateLevels',[levs(2) oslev]);
% width = fin-init
% see attached figure "pulsewidth"

As you can see, these functions can be pretty useful! I believe these functions have the ability to work on trains of pulses too, so you don't even need to epoch your data.

If nothing in here quite works for you, then I am curious what the application of your problem is for. Maybe there is another way to circumvent the issue.

[Eric Kuebler]

I am going to try this in our analysis, and I appreciate the time you 've taken to make the advantages of these functions clear to me.

It is important for our application that the function choses either pre- or post-peak as the lower state consistently. I think that is the only advantage that the code I loaded above may have, I tell it to look in the post-peak portion for that minimum.

The application we are using this for is classification of cell types in single unit recordings of neurons in the visual cortex. The half width duration is one of the parameters that gives good classification.