how to do digital filter on this EEG data

John Williams answered . 2025-03-21 18:40:53

% Given parameters
sample_rate = 50000; % 50 kHz
cutoff = 10; % Cutoff frequency in Hz
filter_order = 500; % Increase filter order for sharper cutoff

% Normalize cutoff frequency
normalized_cutoff = cutoff / (sample_rate / 2);

% Design the low-pass FIR filter
lowpass = fir1(filter_order, normalized_cutoff, 'low');

% Analyze the frequency response of the filter
fvtool(lowpass, 1);

% Apply the filter to your EEG data
y_lowpass = filter(lowpass, 1, y);

% Plot the original and filtered signals
t = (0:length(y)-1) / sample_rate; % Time vector
figure;
subplot(2, 1, 1);
plot(t, y);
title('Original Signal');
xlabel('Time (s)');
ylabel('Amplitude');

subplot(2, 1, 2);
plot(t, y_lowpass);
title('Low-pass Filtered Signal');
xlabel('Time (s)');
ylabel('Amplitude');

Key Points:

1. Filter Order:

   • A higher order (e.g., 500) provides a sharper transition band but increases computation time. You can tune this value based on your needs.

2. Frequency Normalization:

   • Ensure that the cutoff frequency is normalized with respect to the Nyquist frequency ($\text{sample\_rate}/2$).

3. Filter Verification:

   • Use fvtool to visualize the frequency response of the filter and verify that it attenuates frequencies above 10 Hz effectively.

Alternative: Using designfilt

Instead of fir1, you can use designfilt for more control over the filter design:

% Design a low-pass FIR filter with designfilt
d = designfilt('lowpassfir', ...
               'FilterOrder', filter_order, ...
               'CutoffFrequency', cutoff, ...
               'SampleRate', sample_rate);

% Analyze the frequency response
fvtool(d);

% Apply the filter
y_lowpass = filter(d, y);

Verify the Filtered Signal

To ensure that frequencies above 10 Hz are removed:

1. Perform a Fourier Transform of the original and filtered signals to compare the frequency spectra.
2. Use fft and plot to visualize the power in each frequency band.

% FFT of the original signal
Y_orig = fft(y);
f = linspace(0, sample_rate / 2, floor(length(Y_orig) / 2) + 1);

% FFT of the filtered signal
Y_filtered = fft(y_lowpass);

% Plot the magnitude spectra
figure;
plot(f, abs(Y_orig(1:length(f))), 'b');
hold on;
plot(f, abs(Y_filtered(1:length(f))), 'r');
title('Frequency Spectrum');
xlabel('Frequency (Hz)');
ylabel('Magnitude');
legend('Original', 'Filtered');

Troubleshooting

1. If the filtered signal still contains power above 10 Hz, increase the filter order or consider using an IIR filter (e.g., Butterworth) for a sharper cutoff.
2. Ensure that your signal is not being corrupted by noise or artifacts during recording or preprocessing.

Let me know if you need additional guidance!

Not satisfied with the answer ?? ASK NOW