How to remove relative frequency information in time domain?

Hi everyone,

I have an interference signal which contains cross-correlation term and auto-correlation term.

However, there are two signal(46.2 kHz & 69.4 kHz) symmerical to the main signal(57.8 kHz).

Does any method to correct the interference signal for removing/suppressing these two relative frequency(11.6 kHz) in time-domain?

clear all
clc
clf
close all

%%
filename = ('220315_16_27_rawdata.xlsx'); 
data = xlsread(filename);
points = data(:,1);
signal = data(:,2);
signal = detrend(signal(11:1120,1),4);
Fs = 200e3; % sampling rate is 200 kHz
wavelength = linspace(770,830,length(signal));

figure(1)
plot(wavelength,signal,'linewidth',1);
title('commercial spectrum');
xlabel('wavelength (nm)');
ylabel('Amplitude');
set(gca,'fontsize',18);
grid on

load('depth.mat');
window = hanning(length(signal));
signal = window.*signal;
FFT = abs(fftshift(fft(signal,2048)));
FFT = FFT/max(FFT);
logFFT = 10*log(FFT);

figure(2)
plot(depth,logFFT,'linewidth',1);
title('A-scan');
xlabel('Depth (µm)');
ylabel('Amplitude');
xlim([0 3500]);
set(gca,'fontsize',18);
grid on

[Vivian Yu]

I have tried bandpass filter for removing other signal, but it doesn't make sense.

[Star Strider]

I am not certain what that signal represents, however the symmetry suggests double sideband transmitted carrier modulation. See if using the demod function can improve it. (The carrier frequency appears to be 57.8 kHz.)

Answer 1

[Mathieu NOE]

hello

this would be my suggestion, based on using multiple notch filters (simply add a new frequency to the list and the filtering is done automatically).

This will preserve the general aspect of the time signal - a contrario from a bandpass filter that will also remove the static trend

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% select filters to apply
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

option_notch = 1; % 0 = without notch filter , 1 = with notch filter
    fc_notch = [30200; 43300; 45200; 53200; 67400];       % notch center frequencies

option_BPF = 0; % 0 = without bandpass filter , 1 = with bandpass filter
    fc = 56.25e3;      % Hz
    bandwith = 4e3;      % Hz
     N_bpf = 4;           % filter order

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% load signal
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

filename = ('220315_16_27_rawdata.xlsx'); 
data = xlsread(filename);
points = data(:,1);
signal = data(:,2);
% signal = detrend(signal(11:1120,1),4);
Fs = 200e3; % sampling rate is 200 kHz
dt = 1/Fs;
[samples,channels] = size(signal);

% time vector 
time = (0:samples-1)*dt;
signal_filtered = signal;

%% notch filter section %%%%%%
% y(n)=G*[x(n)-2*cos(w0)*x(n-1)+x(n-2)]+[2*p cos(w0)*y(n-1)-p^2 y(n-2)]
if option_notch ~= 0
    for ck =1:numel(fc_notch)
        w0 = 2*pi*fc_notch(ck)/Fs;
        p = 0.925;

        % digital notch (IIR)
        num1z=[1 -2*cos(w0) 1];
        den1z=[1 -2*p*cos(w0) p^2];

        % now let's filter the signal 
        signal_filtered = filter(num1z,den1z,signal_filtered);
    end
end

%% band pass  filter section %%%%%%
if option_BPF ~= 0
    f_low = fc-0.5*bandwith;
    f_high = fc+0.5*bandwith;
    [b,a] = butter(N_bpf,2/Fs*[f_low f_high]);
    signal_filtered = filter(b,a,signal_filtered);
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FFT parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
NFFT = 1000;    % 
OVERLAP = 0.75;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display 1 : time domain plot
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(1),
subplot(211),plot(time,signal,'b');grid on
title(['Time plot  / Fs = ' num2str(Fs) ' Hz / raw data ']);
xlabel('Time (s)');ylabel('Amplitude');
subplot(212),plot(time,signal_filtered,'r');grid on
title(['Time plot  / Fs = ' num2str(Fs) ' Hz / filtered data ']);
xlabel('Time (s)');ylabel('Amplitude');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display 2 : averaged FFT spectrum
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[freq, spectrum_raw] = myfft_peak(signal,Fs,NFFT,OVERLAP);
[freq, spectrum_filtered] = myfft_peak(signal_filtered,Fs,NFFT,OVERLAP);

figure(2),plot(freq,20*log10(spectrum_raw),'b',freq,20*log10(spectrum_filtered),'r');grid on
df = freq(2)-freq(1); % frequency resolution 
title(['Averaged FFT Spectrum  / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(df,3) ' Hz ']);
xlabel('Frequency (Hz)');ylabel('Amplitude (dB (L))');
legend('raw','filtered');

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 percentage of overlap % (between 0 and 0.95)

[samples,channels] = size(signal);

% fill signal with zeros if its length is lower than nfft
if samples<nfft
    s_tmp = zeros(nfft,channels);
    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*ones(1,channels));
        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