Remove harmonics background noise
14 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
Sammy
el 8 de Jul. de 2024
Comentada: Mathieu NOE
el 11 de Jul. de 2024
Hi,
I have a signal sampled at 128 [kHz], it contains noisy spikes at harmonic multiples of 59.6 [Hz]. Tried different bandstop filtering techniques, but can't seem to get a proper result.
Would appreciate it if something with some experience could elaborate how to remove the humming noise in the background.
Cheers
Sammy
2 comentarios
Mathieu NOE
el 8 de Jul. de 2024
hello
this is quite a strange signal - what is it suppose to represent ? looks like a pulse train rather than just a simple corruption by hum noise. How did you come with a 59.6 [Hz] frequency for the spikes ?
looking at different time scales , we have no simple repetitive spikes , these spikes seem to occur more randomly that I expected in first place. Maybe simply a low pass filter would suffice but that depends what you want to keep in your signal
Respuesta aceptada
Mathieu NOE
el 8 de Jul. de 2024
I have some doubts we can recover a clean sinewave from that ... nevertheless I tried a few things ... first I decimated the signal by a factor of 20 as Fs = 128kHz is pretty fast for a signal that does not exceed 2 kHz
then bandpass filtered it - the spectrogram tells us that there is indeed some kind of sweeip burried in the noise , but it starts above 200 Hz obviously and also its amplitude is not constant with time
here we are for the moment :
then we will try to remove the horizontal lines with a notch filter - looped as many times as there are H lines (harmonics) .
But we need to make sure the frame rate is measures as accurately as possible. So I picked the peaks of the 1 second of signal and theit time indexes are checked to be regurarly spaced . But the computed frequency here is 62.0456 hz and not 59.6 Hz (this is due to the lack of frequency resolution of your fft computation). So with that info , we can run a notch filter a this frequency and its N harmonics
the new spectrogram appears "cleaned" from these H lines , but there is still an important background noise, so the time signal is still very noisy.
NB : my code uses the function peakseek PeakSeek - File Exchange - MATLAB Central (mathworks.com) that I provide in attachment. A simpler yet very effective alternative to findpeaks !
clc
clearvars
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% load signal
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
load('short_sweep.mat')
signal = s(:); % 1 channel of data
clear s
dt = 1/128e3; % sampling time interval
Fs = 1/dt; % sampling rate is 128kHz here
[samples,channels] = size(signal);
t = (0:samples-1)*dt;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% compute "clicks" period on first second
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
t1 = t(1:Fs);
signal1 = signal(1:Fs);
% try with peakseek
minpeakdist = Fs/100; % in samples
minpeakh = 2;
[locs1, pks1]=peakseek(signal1,minpeakdist,minpeakh);
% plot(t1,signal1,t1(locs1),signal1(locs1),'dk');
period = mean(diff(t1(locs1)));
ff = 1/period;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% decimate (if needed)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% NB : decim = 1 will do nothing (output = input)
decim = 20;
if decim>1
Fs = Fs/decim;
for ck = 1:channels
newsignal(:,ck) = decimate(signal(:,ck),decim);
end
signal = newsignal;
end
samples = length(signal);
time = (0:samples-1)*1/Fs;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% band pass filter section %%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
f_low = 100;
f_high = 2500;
N_bpf = 2; % filter order
[b,a] = butter(N_bpf,2/Fs*[f_low f_high]);
signal = filtfilt(b,a,signal);
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% notch filter
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fc_notch = ff; % 62.0456 Hz in fact
rho = 0.997; % something slightly less than 1 (will change the width and the depth of the notch filter).
N = 45; % how many harmonics of fc_notch do we want to remove ?
for k = 1:N+1
w0 = 2*pi*k*fc_notch/Fs;
% digital notch (IIR)
num1z=[1 -2*cos(w0) 1];
den1z=[1 -2*rho*cos(w0) rho^2];
% now let's filter the signal
signal = filtfilt(num1z,den1z,signal);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% FFT analysis
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
nfft = 512;
overlap = 0.75;
dB_range = 40; % display this dB range on y axis
[S,F,T] = myspecgram(signal, Fs, nfft, overlap); % overlap is 75% of nfft here
sg_dB = 20*log10(abs(S)); % expressed now in dB
% saturation of the dB range : the lowest displayed level is spectrogram_dB_scale dB below the max level
min_disp_dB = round(max(max(sg_dB))) - dB_range;
%% plots
figure(1);
plot(time,signal)
figure(2);
imagesc(T,F,sg_dB)
caxis([min_disp_dB min_disp_dB+dB_range]);
% ylim([100 Fs/2.56]);
hcb = colorbar('vert');
set(get(hcb,'Title'),'String','dB')
set(gca,'YDir','Normal')
xlabel('Time (secs)')
ylabel('Freq (Hz)')
title('Specgram')
colormap('jet');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [ZxP,ZxN] = find_zc(x,y,threshold)
% put data in rows
x = x(:);
y = y(:);
% positive slope "zero" crossing detection, using linear interpolation
y = y - threshold;
zci = @(data) find(diff(sign(data))>0); %define function: returns indices of +ZCs
ix=zci(y); %find indices of + zero crossings of x
ZeroX = @(x0,y0,x1,y1) x0 - (y0.*(x0 - x1))./(y0 - y1); % Interpolated x value for Zero-Crossing
ZxP = ZeroX(x(ix),y(ix),x(ix+1),y(ix+1));
% negative slope "zero" crossing detection, using linear interpolation
zci = @(data) find(diff(sign(data))<0); %define function: returns indices of +ZCs
ix=zci(y); %find indices of + zero crossings of x
ZeroX = @(x0,y0,x1,y1) x0 - (y0.*(x0 - x1))./(y0 - y1); % Interpolated x value for Zero-Crossing
ZxN = ZeroX(x(ix),y(ix),x(ix+1),y(ix+1));
end
function [fft_specgram,freq_vector,time] = myspecgram(signal, Fs, nfft, Overlap)
% FFT peak spectrogram of signal (example sinus amplitude 1 = 0 dB after fft).
% signal - input signal,
% Fs - Sampling frequency (Hz).
% nfft - FFT window size
% Overlap - buffer overlap % (between 0 and 0.95)
signal = signal(:);
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_specgram = [];
for ci=1:spectnum
start = (ci-1)*offset;
sw = signal((1+start):(start+nfft)).*window;
fft_specgram = [fft_specgram abs(fft(sw))*4/nfft]; % X=fft(x.*hanning(N))*4/N; % hanning only
end
% 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_specgram = fft_specgram(select,:);
freq_vector = (select - 1)*Fs/nfft;
% time vector
% time stamps are defined in the middle of the buffer
time = ((0:spectnum-1)*offset + round(nfft/2))/Fs;
end
4 comentarios
Mathieu NOE
el 9 de Jul. de 2024
my pleasure !
if my answer has helped you , maybe you could consider accepting it ! :)
Mathieu NOE
el 11 de Jul. de 2024
fyi , some digital filters matlab code in attachment
Más respuestas (0)
Ver también
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!