Remove harmonics background noise

Question

Sammy 2024-7-8，7:41

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2135376-remove-harmonics-background-noise

评论： Mathieu NOE 2024-7-11，8:43

short_sweep.mat

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 个评论
显示无隐藏无

Mathieu NOE 2024-7-8，8:48

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

Sammy 2024-7-8，9:00

编辑：Sammy 2024-7-8，9:05

This signal was acquired using a rolling shutter camera which has a very short row time of 8 [us] and a framerate of 59.6 FPS, which is the source of the noise. The signal itself contains the humming noise and in the background a sine sweep signal which should range between 200-2000 [Hz].
The spikes I was referring to are in the frequency domain,(allthough there are also some spikes in the time domain), which you can see if you look at its PSD:

figure;

S=fft(short_sweep);

power_spectrum = abs(S).^2;

frequencies = linspace(0, Fs/2, length(power_spectrum)/2 + 1);

plot(frequencies,power_spectrum(1:length(frequencies)));

title('Power Spectrum of s(t) (resampled) - unfiltered');

xlim([0 10000]);

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Mathieu NOE 2024-7-8，13:57

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2135376-remove-harmonics-background-noise#answer_1482811

在 MATLAB Online 中打开

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 个评论
显示 2更早的评论隐藏 2更早的评论

Mathieu NOE 2024-7-8，15:56

在 MATLAB Online 中打开

here some improvement;

I added a time varying bandpass filter , now the signal is closer to a chirp only (sure with some amplitude variations and a bit of noise)

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);
dt = 1/Fs;
time = (0:samples-1)*dt;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% time varying IIR band pass filtering 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
BW = 100;
fc_init = BW+50;
slop = 2150/5.6;  % Hz/s
y = signal; % init y with input signal1
N = 3; % filter's order
nB = 2*N+1;
nA = 2*N+1;
for k = nB:samples % for iteration >= nB sample
    fc = max(fc_init,slop*k*dt);
    flowk = fc-BW/2;
    fhighk = fc+BW/2;
    [B,A] = butter(N,2/Fs*[flowk fhighk]);
    % main recursion
    y(k) = B*signal(k:-1:k-nB+1) - A(2:nA)*y(k-1:-1:k-nA+1);
end
signal(0.25*Fs:samples) = y(0.25*Fs:samples); % the beginning of the filtered signal show an important ringing due to the butterworth transient 
                                               % so a dumb solution is to replace the first filtered samples by the input signal
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% FFT analysis
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
nfft = 512; 
overlap = 0.75; 
dB_range = 40; % display this dB range on Z 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