2nd order butterworth filter and PCA

Good morning, I am performing a FFT on my accelerometer data, but I am having issues with the 2nd-order Butterworth filter and pca interpretation. Please find a snippet of my code below.

RE the 2nd-order Butterworth filter, I am trying to first, downsample to 100 Hz (haven't managed to do this yet), and then bandpass filter using a second-order Butterworth filter between 3-15 Hz. However, I keep getting the error: "The cutoff frequencies must be within the interval of (0,1). Would anyone be able to assist me with this - specifically, highlight where in the code I am going wrong?

Regarding the pca, the anwser below (bolded) is the anwser from the pca. The columns represent 'x, z, y' - does this suggest that the Y is the most dominant (movement-centred) axis of tremor acceleration? From plotting, this seems correct but I'd like some confirmation (plot not included in the code provided below).

PCA ans = 0.9780 -0.1221 0.1694

-0.0357 0.7013 0.7120

0.2057 0.7023 -0.6815

Code:

t = data(:,1); % column corresponding to time (block of 60s)

Accel = data(:,2:end); % three columns corresponding to x, y, z (accelerometer data)

srate = 5000; % Hz (need to figure out how to downsample to 100Hz)

Fs = 1/srate; % Sampling Rate (in Hz)

Fn =Fs/2; % Nyquist Frequency (in Hz)

coeff = pca(Accel); % Need to obtain the most dominant (movement-centred) axis of tremor acceleration; use first component for all subsequent analyses

%2nd order Butterworth filter between 3-15 Hz

Wp = [3 15]/Fn; % Norrmalised Passband Frequencies

Ws = [1.5 15.99]/Fn; % Normalised Stopband Frequencies

Rp = 1; % Passband Ripple (dB)

Rs = 25; % Stopband Ripple (dB)

[n,Wn] = buttord(Wp, Ws, Rp, Rs); % Optimal Filter Order

[b,a] = butter(n, Wn); % Calculate Filter Coefficients

[sos,g] = tf2sos(b,a); % Convert To Second-Order Sections For Stability

figure(1)

freqz(sos, 4096, Fs) % Filter Bode Plot

Many thanks in advance for any assistance :)

Answer 1

[Star Strider]

With respect to downsampling, use the resample function. I wouldn’t downsample the signal unless you absolutely must for some reason. Most of the time, I use resample to create a unifformly-sampled signal from a non-uniformly-sampled signal, since all signal processing requires uniform sampling. If you do decide to downsample, be sure to change ‘Fs’ to reflect that, if you’re doing the filtering after the downsampling.

I don’t understand this assignmment:

srate = 5000;

If ‘srate’ is your sampling rate, it’s actually ‘Fs’, the sampling frequency, so assign it to ‘Fs’ instead. That’s also the indirect cause of the error you’re getting.

With that correction, the rest of your code should work. (It looks like I wrote it!)

I would also make two change these assignments:

[b,a] = butter(n, Wn);                                   % Calculate Filter Coefficients
[sos,g] = tf2sos(b,a);                                   % Convert To Second-Order Sections For Stability

to:

[z.p,k]] = butter(n, Wn);                                % Calculate Filter Coefficients
[sos,g] = zp2sos(b,a);                                   % Convert To Second-Order Sections For Stability

That will design a better filter that the code I originally wrote.

And of course, use filtfilt to do the actual filtering.

[KR]

Thank you for such a prompt response! I am glad you responded, because this is your code! I have created most of my code based on your assistance from other questions on the forum (so thank you for providing such understandable code!).

RE the downsampling, I don't think it is necessary, but I was basing it on what is outlined in previous literature investigating a similar research question to mine. RE strate, I have replaced this to 'Fs' - thanks for this tip :)

I have applied the recommended changes (below) to replace the above mentionned code, but keep getting 'undefined function or variable 'b'. I have tried to replace 'b, a' with 'p,k' and that doesn't seem to work.

Any thoughts as to why this is occuring? Also, regarding the filtfilt, do you mean to implent a line of code after the butterworth, as I did in the code below?

Thanks again for your help - I greatly appreciate it :)!

[z.p,k] = butter(n, Wn);                                % Calculate Filter Coefficients
[sos,g] = zp2sos(b,a);                                   % Convert To Second-Order Sections
y = filtfilt(sos,g, Accel);

[Star Strider]

As always, my pleasure!

I very much appreciate your compliment!

If there is no specific need to downsample, I would reccommend against doing it. I don’t know what you’re doing, however more data is always better.

That was a typo in the conversion to the zero-pole-gain representation (it’s been a long day).

That line should be:

[sos,g] = zp2sos(z,p,k);                                % Convert To Second-Order Sections For Stability

The code should now work correctly.

[KR]

Yes, it's been lovely (indirectly) learning from you! Thanks again :)

Success - it is working!! Sorry - I should have caught that typo (that's my fault)!

Also, if you have the opportunity, do you mind highlighting the error in this code:

[pxx,f] = pwelch(Accel,hanning(300000),[],300000,5000);  
[pxx,f] = periodogram(Accel,300000,5000,hanning(300000));

I am trying to calculate the peak frequency between 3-15 Hz using a sliding window of 1 s over each 3-s window with 50% overlap (the block is 60 s). Above is the code I am running with functions 'pwelch' and 'periodogram' (prefer to use periodogram but it does not seem to be working!). The periodogram error is suggesting that the sampling frequency should be a scalar (which I am pretty sure it is). Any thoughts? Worth noting, that I believe that '5000' (the sample frequency) should be replaced with the nyquist (i.e. 2500). Thanks again for your help :)