MATLAB: How to make Narrow Bandpass filter

bandpassfilter

I have a signal, that has a bunch of noise and sine signals in it, 4 sine signals too be exact, and 4 noise signals.

i have pin pointed the frequencies of these sine signals to be 159.0574, 359.3519, 459.4992, 709.8674.

The signal and its spektrum looks like this:

I need to filter out all the noise, so i tried too use bandpass, using each sine frequency :

f1 = bandpass(signalas, [floor(virsunes(1)) ceil(virsunes(1))], Fd);

f2 = bandpass(signalas, [floor(virsunes(2)) ceil(virsunes(2))], Fd);

f3 = bandpass(signalas, [floor(virsunes(3)) ceil(virsunes(3))], Fd);

f4 = bandpass(signalas, [floor(virsunes(4)) ceil(virsunes(4))], Fd);

virsunes is just an array of those frequencies.

When i combined the filtered signals, i got this:

A lot of the noise has stuck around because it's too close too the sine signal. How to make a bandpass filter more narrow, for it too filter all the noise?

I attached the signal file "G4_12.dot".

Best Answer

hello

seems I have already seen similar post somewhere... a buddy of yours ?

you did not mentionned the sampling frequency , but 24 kHz seemed correct

here my little contribution :


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%





% options 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% spectrogram dB scale
spectrogram_dB_scale = 60;  % dB range scale (means , the lowest displayed level is XX dB below the max level)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% load signal
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% test data
Fs = 24000;
% t = 0:1/Fs:10-1/Fs;
fid = fopen('G4_12.dat','r')
[signal,samples] = fscanf(fid, '%12f');
fclose(fid);
% decimate
decim = 12;
signal = decimate(signal,decim);
Fs = Fs/decim;
samples = length(signal);
time =(1:samples)/Fs;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FFT parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
n_max = fix(log2(samples));
NFFT = 2^(n_max-1);    % 2^n with n chosen so that NFFT is below length of signal (even after decimation) by factor > 2 so some averaging is doable
% NFFT = 2048/4;    % 2^n with n chosen so that NFFT is below length of signal (even after decimation) by factor > 2 so some averaging is doable
NOVERLAP = round(0.95*NFFT);
w = hanning(NFFT);     % Hanning window / Use the HANN function to get a Hanning window which has the first and last zero-weighted samples.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



% display : averaged FFT spectrum before filter
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[sensor_spectrum, freq] = pwelch(signal,w,NOVERLAP,NFFT,Fs,'power');
sensor_spectrum_dB = 10*log10(sensor_spectrum);% convert to dB scale (ref = 1)

% sensor_spectrum_dB = sensor_spectrum_dB-min(sensor_spectrum_dB); % y axis shift so dB values are positive (min value = 0)
%% findpeaks
% df = Fs/NFFT;
% MAXW = 10*df;
% MPD = 20;
% MPP = 20;
% [peaks,loc,W,P] = findpeaks(sensor_spectrum_dB,'MinPeakDistance',MPD,'MaxPeakWidth',MAXW,'MinPeakProminence',MPP);
% freq_peaks = freq(loc);
freq_peaks = [159.0574, 359.3519, 459.4992, 709.8674];
peaks = interp1(freq,sensor_spectrum_dB,freq_peaks,'linear');
figure(1),plot(freq,sensor_spectrum_dB,'b',freq_peaks,peaks,'*r');grid
title(['Averaged FFT Spectrum  / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(freq(2)-freq(1)) ' Hz ']);
xlabel('Frequency (Hz)');ylabel(' dB')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% display 2 : time / frequency analysis : spectrogram demo
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[sg,fsg,tsg] = specgram(signal,NFFT,Fs,hanning(NFFT),NOVERLAP);  
% FFT normalisation and conversion amplitude from linear to dB (peak)
sg_dBpeak = 20*log10(abs(sg))+20*log10(2/length(fsg));     % NB : X=fft(x.*hanning(N))*4/N; % hanning only
% scalingf the dB range : 
min_disp_dB = round(max(max(sg_dBpeak))) - spectrogram_dB_scale;
sg_dBpeak(sg_dBpeak<min_disp_dB) = min_disp_dB;
% plots spectrogram
figure(2);
imagesc(tsg,fsg,sg_dBpeak);colormap('jet');
axis('xy');colorbar('vert');grid
title(['Spectrogram / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(fsg(2)-fsg(1)) ' Hz ']);
xlabel('Time (s)');ylabel('Frequency (Hz)');
% keep only signal content extracted by bandpass filters 
fc = freq_peaks;      % Hz

bandwith = 8;      % Hz
f_low = fc-0.5*bandwith;
f_high = fc+0.5*bandwith;
N = 4;
signal_filtered = zeros(size(signal));
for ci = 1:length(fc)
    [b,a] = butter(N,2/Fs*[f_low(ci) f_high(ci)]);
    % add filtered signal to total final signal
    tmp = filtfilt(b,a,signal);
    signal_filtered = signal_filtered + tmp;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display : averaged FFT spectrum before / after notch filter
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[sensor_spectrum_filtered, freq] = pwelch(signal_filtered,w,NOVERLAP,NFFT,Fs,'power');
sensor_spectrum_filtered_dB = 10*log10(sensor_spectrum_filtered);% convert to dB scale (ref = 1)
figure(4),plot(freq,sensor_spectrum_dB,'b',freq,sensor_spectrum_filtered_dB,'r');grid
title(['Averaged FFT Spectrum  / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(freq(2)-freq(1)) ' Hz ']);
legend('before filter','after filter');
xlabel('Frequency (Hz)');ylabel(' dB')
figure(5),plot(time,signal_filtered,'b');grid
title('Filtered signal');
xlabel('Time (s)');ylabel(' amplitude')

Related Solutions

MATLAB: How to only pass the signals that I want

hello

this code will do fft / spectrogram and bandpass filtering - adapt the corner frequencies to your own needs !

clc
clear all
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%





% FFT parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
NFFT = 4096;    % 
OVERLAP = 0.5;
% spectrogram dB scale
spectrogram_dB_scale = 80;  % dB range scale (means , the lowest displayed level is XX dB below the max level)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% options 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% if you are dealing with acoustics, you may wish to have A weighted
% spectrums 
% option_w = 0 : linear spectrum (no weighting dB (L) )
% option_w = 1 : A weighted spectrum (dB (A) )
option_w = 1;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% load signal
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% data
% [data,Fs] = audioread('myWAVaudiofile.wav');
[data,Fs] = audioread('Immigran_Hiss.wav');
channel = 1;
signal = data(:,channel);
samples = length(signal);
%% decimate (if needed)
% NB : decim = 1 will do nothing (output = input)
decim = 1;
if decim>1
    signal = decimate(signal,decim);
    Fs = Fs/decim;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% display 1 : averaged FFT spectrum
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[freq, sensor_spectrum] = myfft_peak(signal,Fs,NFFT,OVERLAP);
% convert to dB scale (ref = 1)
sensor_spectrum_dB = 20*log10(sensor_spectrum);
% apply A weigthing if needed


if option_w == 1
    pondA_dB = pondA_function(freq);
    sensor_spectrum_dB = sensor_spectrum_dB+pondA_dB;
    my_ylabel = ('Amplitude (dB (A))');
else
    my_ylabel = ('Amplitude (dB (L))');
end
figure(1),plot(freq,sensor_spectrum_dB,'b');grid
title(['Averaged FFT Spectrum  / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(freq(2)-freq(1)) ' Hz ']);
xlabel('Frequency (Hz)');ylabel(my_ylabel);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



% display 2 : time / frequency analysis : spectrogram demo
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[sg,fsg,tsg] = specgram(signal,NFFT,Fs,hanning(NFFT),floor(NFFT*OVERLAP));  
% FFT normalisation and conversion amplitude from linear to dB (peak)

sg_dBpeak = 20*log10(abs(sg))+20*log10(2/length(fsg));     % NB : X=fft(x.*hanning(N))*4/N; % hanning only

% apply A weigthing if needed
if option_w == 1
    pondA_dB = pondA_function(fsg);
    sg_dBpeak = sg_dBpeak+(pondA_dB*ones(1,size(sg_dBpeak,2)));
    my_title = ('Spectrogram (dB (A))');
else
    my_title = ('Spectrogram (dB (L))');
end
% saturation of the dB range : 

% saturation_dB = 60;  % dB range scale (means , the lowest displayed level is XX dB below the max level)

min_disp_dB = round(max(max(sg_dBpeak))) - spectrogram_dB_scale;
sg_dBpeak(sg_dBpeak<min_disp_dB) = min_disp_dB;
% plots spectrogram

figure(2);
imagesc(tsg,fsg,sg_dBpeak);colormap('jet');
axis('xy');colorbar('vert');grid
title([my_title ' / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(fsg(2)-fsg(1)) ' Hz ']);
xlabel('Time (s)');ylabel('Frequency (Hz)');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% removal of hiss between t = 0 and t = 1.68 s
% use of bandpass filter to remove noise below 500 and above 2500 Hz
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
data_filtered = data;
t_min = 0;
t_max = 1.68;
ind = 1+fix(t_min*Fs:t_max*Fs-1);
% keep only signal content extracted by bandpass filters 
f_low = 500;
f_high = 2500;
N = 4;
% signal_filtered = zeros(size(signal));
[b,a] = butter(N,2/Fs*[f_low f_high]);
% now let's filter the signal
data_filtered(ind,:) = filtfilt(b,a,data(ind,:));
% save as wav file 
audiowrite('Immigran_Hiss-filt.wav',data_filtered,Fs);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display 2B : time / frequency analysis : spectrogram demo
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[sg,fsg,tsg] = specgram(data_filtered(:,1),NFFT,Fs,hanning(NFFT),floor(NFFT*OVERLAP));  
% FFT normalisation and conversion amplitude from linear to dB (peak)
sg_dBpeak = 20*log10(abs(sg))+20*log10(2/length(fsg));     % NB : X=fft(x.*hanning(N))*4/N; % hanning only
% apply A weigthing if needed
if option_w == 1
    pondA_dB = pondA_function(fsg);
    sg_dBpeak = sg_dBpeak+(pondA_dB*ones(1,size(sg_dBpeak,2)));
    my_title = ('Spectrogram (dB (A))');
else
    my_title = ('Spectrogram (dB (L))');
end
% saturation of the dB range : 
% saturation_dB = 60;  % dB range scale (means , the lowest displayed level is XX dB below the max level)
min_disp_dB = round(max(max(sg_dBpeak))) - spectrogram_dB_scale;
sg_dBpeak(sg_dBpeak<min_disp_dB) = min_disp_dB;
% plots spectrogram
figure(3);
imagesc(tsg,fsg,sg_dBpeak);colormap('jet');
axis('xy');colorbar('vert');grid
title([my_title ' / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(fsg(2)-fsg(1)) ' Hz ']);
xlabel('Time (s)');ylabel('Frequency (Hz)');
function pondA_dB = pondA_function(f)
	% dB (A) weighting curve
	n = ((12200^2*f.^4)./((f.^2+20.6^2).*(f.^2+12200^2).*sqrt(f.^2+107.7^2).*sqrt(f.^2+737.9^2)));
	r = ((12200^2*1000.^4)./((1000.^2+20.6^2).*(1000.^2+12200^2).*sqrt(1000.^2+107.7^2).*sqrt(1000.^2+737.9^2))) * ones(size(f));
	pondA = n./r;
	pondA_dB = 20*log10(pondA(:));
end
function  [freq_vector,fft_spectrum] = myfft_peak(signal, Fs, nfft, Overlap)
% FFT peak spectrum of signal  (example sinus amplitude 1   = 0 dB after fft).
% Linear averaging
%   signal - input signal, 
%   Fs - Sampling frequency (Hz).
%   nfft - FFT window size
%   Overlap - buffer overlap % (between 0 and 0.95)
samples = length(signal);
% fill signal with zeros if its length is lower than nfft
if samples<nfft
    s_tmp = zeros(nfft,1);
    s_tmp((1:samples)) = signal;
    signal = s_tmp;
    samples = nfft;
end
% window : hanning
window = hanning(nfft);
window = window(:);
%    compute fft with overlap 
 offset = fix((1-Overlap)*nfft);
 spectnum = 1+ fix((samples-nfft)/offset); % Number of windows
%     % for info is equivalent to : 
%     noverlap = Overlap*nfft;
%     spectnum = fix((samples-noverlap)/(nfft-noverlap));	% Number of windows
    % main loop
    fft_spectrum = 0;
    for i=1:spectnum
        start = (i-1)*offset;
        sw = signal((1+start):(start+nfft)).*window;
        fft_spectrum = fft_spectrum + (abs(fft(sw))*4/nfft);     % X=fft(x.*hanning(N))*4/N; % hanning only 
    end
    fft_spectrum = fft_spectrum/spectnum; % to do linear averaging scaling
% one sidded fft spectrum  % Select first half 
    if rem(nfft,2)    % nfft odd
        select = (1:(nfft+1)/2)';
    else
        select = (1:nfft/2+1)';
    end
fft_spectrum = fft_spectrum(select);
freq_vector = (select - 1)*Fs/nfft;
end

MATLAB: Efficient pass filters for EEG signals

The lowpass function (and those like it) have 2 outputs, the second being a digital filter object. Save the digital filter object and use it with filtfilt to use it with other signals.

That is much more efficient than designing it each time you use it, which is what it appears you are doing now. So, design it once and use it as often as necessary afterwards.

Best Answer

Related Solutions

MATLAB: How to only pass the signals that I want

MATLAB: Efficient pass filters for EEG signals

Related Question