Need help designing stopband filter

Hello,

I'm trying to design a stopband filter that can remove spikes at 1.294 Hz and 3.101 Hz

My data matrix 64x64x64x2960, where the three first dimensions are spatial dimensions and the 4th dimension is the temporal dimension along which I need to filter the data.

I have tried this, but it doesnt work for 2960 data points and requires at least 5040 data points:

Fs = 10;

fcomb = [[1.285 1.29 1.3 1.35],[2.85 2.90 3.15 3.20]];

mags = [[1 0 1], [0 1];

dev = [[0.5 0.1 0.5], [0.1 0.5]];

[n,Wn,beta,ftype] = kaiserord(fcomb,mags,dev,Fs);

hh = fir1(n,Wn,ftype,kaiser(n+1,beta),'noscale');

figure

freqz(hh, 1, 2^20, Fs)

set(subplot(2,1,1), 'XLim',[45 65]) % Optional

set(subplot(2,1,2), 'XLim',[45 65])

I have also tried matlabs bandstop() function, but somehow it made my data completely flat in freq domain. I only want to remove the spikes at 1.294 Hz and 3.101 Hz ...

Answer 1

[Star Strider]

Thank you for quoting my code!

The FIR filter will eliminate both frequencies in one filtfilt call to it, however they are not efficient and can only be used on long signals. It might be best to use two IIR filters in series.

That design would be something like this —

Fs = 10;
Fn = Fs/2;
Ws = [1.290   1.298]/Fn;
Wp = Ws.*[0.999  1.001];
Rs = 50;
Rp =  1;
[n,Wn] = ellipord(Wp,Ws,Rp,Rs);
[z,p,k] = ellip(n,Rp,Rs,Wp,'stop');
[sos,g] = zp2sos(z,p,k);

figure
freqz(sos, 2^16, Fs)
set(subplot(2,1,1), 'XLim',[1.25 1.35])
set(subplot(2,1,2), 'XLim',[1.25 1.35])

Ws = [3.0 3.2]/Fn;
Wp = Ws.*[0.999 1.001];
Rs = 50;
Rp =  1;
[n,Wn] = ellipord(Wp,Ws,Rp,Rs);
[z,p,k] = ellip(n,Rp,Rs,Wp,'stop');
[sos,g] = zp2sos(z,p,k);

figure
freqz(sos, 2^16, Fs)
set(subplot(2,1,1), 'XLim',[2.5 3.5])
set(subplot(2,1,2), 'XLim',[2.5 3.5])

There is no way of cascading them functionally, so the only option is to use the output of the first filter as the input to the second filter. Use the filtfilt function to use them to filter the signals.

.

[Pseudoscientist]

Thank you very much!

[Star Strider]

As always, my pleasure!

[Pseudoscientist]

Curisouly, even after applying the filter, the spike remains in the data hmm..

My data is a 4D matrix, 64x64x64x2960, where 2960 is the temporal dimension along which i need to filter, so I apply the filter like this:

y = filtfilt(sos,g,permute(double(timeseries),[4 1 2 3]));,

where timeseries is the 64x64x64x2960 matrix that I permute to a 2960x64x64x64 matrix

[Star Strider]

From the documentation for x-input signal :

• filtfilt operates along the first array dimension of x with size greater than 1.

For a 2-D matrix, this means that it operates across dimension 1 (rows), so down each column, and N-D arrays are permitted. With a 2-D matrix,each row must be assigned a different time instant, and each column represents a different signal sampled at the same time instants.

I am not certain which of the two filters is being used here.

Since I do not understand how the matrix is constructed, my only suggestion is to see if it works with the original matrix first. Then filter across one dimension at a time until it gives the desired result, and then if necessary, use permute to give the desired result for this and any similar matrices.

I will help as I can, however I remain a bit mystfied as to what is being filtered.

Also, just to be clear, I should have made the following change to the original filter code:

[sos{1},g{1}] = zp2sos(z,p,k);                                  % First Filter

and:

[sos{2},g{2}] = zp2sos(z,p,k);                                  % Second Filter

Otherwise, calculating them in the same code as in the code I posted, will over-write the first filter coefficient matrix and gain with the second filter coefficient matrix and gain, so only the second frequency (3.101 Hz) will be filtered. (They don’t have to be stored as cell arrays, however this makes addressing them easier. They just need to be distinguished from each other.) I am not certain how your code is written, so this may not be a problem. The filters cannot be functionally cascaded, so the only option is to filter the output of the first filter with the second filter to do the complete filtering, and the filter coefficient matrices and gains need to be distinguished in order for that to work correctly.

.

[Pseudoscientist]

I will try this.

The in the 4D matrix, 64x64x64 are spatial dimensions, so you can imagine there are 2960 three-dimensional cubes. 2960 timepoints of the cube. I need to filter the values of the cube in temporal dimension, so there are 64*64*64 = 262144 signals to filter. For example, this is one signal that I need to filter: timeseries(1,1,1,:), and so on for every x, y, and z from 1 to 64.

[Pseudoscientist]

I have attached an example signal retains the 1.294 Hz spike after filtering it.

[Pseudoscientist]

Heres the code

load('example_signal.mat');

Fs = 10;

Fn = Fs/2;

Ws = [1.290 1.298]/Fn;

Wp = Ws.*[0.999 1.001];

Rs = 50;

Rp = 1;

[n,Wn] = ellipord(Wp,Ws,Rp,Rs);

[z,p,k] = ellip(n,Rp,Rs,Wp,'stop');

[sos1,g1] = zp2sos(z,p,k);

vida_mean_signal_xy_filt = filtfilt(sos1,g1,vida_mean_signal_xy);

N1_xy=2^nextpow2(length(vida_mean_signal_xy_filt));

xy_residual_fft=fft(vida_mean_signal_xy_filt,N1_xy);

xy_residual_fft=xy_residual_fft(1:length(xy_residual_fft)/2);%Discard Half of Points

xy_residual_fft_mag=abs(xy_residual_fft); %Magnitude Spectrum

xy_residual_fft_phase=angle(xy_residual_fft); %Phase Spectrum

fs=10;

f=fs*(0:length(xy_residual_fft)-1)/length(xy_residual_fft); %Frequency axis

figure(1);hold on;clf

plot(f,(xy_residual_fft_mag)); %Plotting the Magnitude Spectrum after Normalization

xlabel('Frequency (Hz)');

ylabel('Magnitude Spectrum');

title('Spectrum')

xlim([0.01 3.4])

[Star Strider]

I cannot run the code, however I see no problems with it.

One thing to note is that the length function returns the size of the largest dimension. The size funciton allows better control over that result.

I will look at the matrix and filtering in a few minutes.

[Pseudoscientist]

Good point, this case the example is 1D luckily

[Star Strider]

I just recently saw the .mat file. My apologies for the resulting delay.

I beleive the frequencies are not defined correctly. The sampling frequency implies the existence of an associated time vector, and so I created one. That information is required in order to design the filter correctly.

This is the most recently posted code, with my corrections, followed by a plot figure of the signal spectrum before and after filtering —

LD = load('example_signal.mat');
vida_mean_signal_xy = LD.vida_mean_signal_xy;                                       % Signal Vector
Fs = 10;
L = numel(vida_mean_signal_xy);                                                     % Signal Vector Length
t = linspace(0, L-1, L)/Fs;                                                         % Time Vector

Fs = 10;
Fn = Fs/2;
Ws = ([-0.01  0.01]+0.646973)/Fn;
Wp = Ws.*[0.99  1.01];
Rs = 50;
Rp =  1;
[n,Wn] = ellipord(Wp,Ws,Rp,Rs);
[z,p,k] = ellip(n,Rp,Rs,Wp,'stop');
[sos1,g1] = zp2sos(z,p,k);

figure
freqz(sos1, 2^18, Fs)
set(subplot(2,1,1), 'XLim',[0 1.3])
set(subplot(2,1,2), 'XLim',[0 1.3])

vida_mean_signal_xy_filt = filtfilt(sos1,g1,vida_mean_signal_xy);
N1_xy=2^nextpow2(L);
xy_residual_fft=fft([vida_mean_signal_xy_filt; vida_mean_signal_xy].',N1_xy)/L;
Fv = linspace(0, 1, N1_xy/2+1)*Fn;                                                  % Frequency Vector
Iv = 1:numel(Fv);                                                                   % Index Vector
% xy_residual_fft=xy_residual_fft(1:length(xy_residual_fft)/2);%Discard Half of Points
% xy_residual_fft_mag=abs(xy_residual_fft); %Magnitude Spectrum
% xy_residual_fft_phase=angle(xy_residual_fft); %Phase Spectrum
% fs=10;
% f=fs*(0:length(xy_residual_fft)-1)/length(xy_residual_fft); %Frequency axis

figure
plot(t, vida_mean_signal_xy)
hold on
plot(t, vida_mean_signal_xy_filt)
hold off
% plot(f,(xy_residual_fft_mag)); %Plotting the Magnitude Spectrum after Normalization
xlabel('Frequency (Hz)');
ylabel('Magnitude Spectrum');
title('Spectrum')
% xlim([0.01 3.4])

figure
subplot(2,1,1)
plot(Fv, abs(xy_residual_fft(Iv,2)*2))
grid
title('Original Signal Spectrum')
subplot(2,1,2)
plot(Fv, abs(xy_residual_fft(Iv,1)*2))
grid
xlabel('Frequency (Hz)')
title('Filtered Signal Spectrum')

Sos the 3.101 Hz filter appears to have worked correctly. The 1.294 Hz appears not to have been specified correctly. (Either that, or it worked correctly, because I see nothing at that frequency.) Specifying it instead as 0.647 Hz works, assuming we are discussing the same frequency spike in the Fourier transform.

.

[Pseudoscientist]

Thank you, it appears my frequency vector in incorrect, sorry for the delayed reply, I was having corona

[Star Strider]

As always, my pleasure!

No worries!

I had it too in March 2020. Not fun, although no serious problems. I hope you’re O.K. now, with no lingering effects.

[Pseudoscientist]

The virus left my cognitive abilities slightly lowered; this may result in me having to post more questions here. Haha.

[Star Strider]

Noi worries about the questions!

Otherwise, I well know that feeling! SARS-CoV-2 precipitates a sustained immune response with elevated cytokines and interleukins of several varieties. That’s likely the cause of the ‘brain fog’, and it takes a while to get back to normal.