How to seperate a superimposed sinusoidal wave from a signal

Could you assist me in isolating a sinusoidal wave that is overlaid(superimposed) on a signal?

I've attached a .mat file where a sine wave is superimposed between the 613th and 1569th data points, which have a sampling frequency of 40Hz. My goal is to extract this sine wave from the entire signal. I've attempted 'blind source separation' techniques, FFT, but it won't works. the reason might be the frequency of introduced sine wave is too small,only 1Hz.

Additionally, could you provide suggestions on how to identify the point at which a sinusoidal wave is introduced if its location is unknown? Your help would be greatly appreciated.

Thank you in advance!

Comment if any other informantion needed.

%%

Answer 1

[Mathieu NOE]

hello

try the code provided below

I found that the sinus signal frequency is probably around 0.5 and not 1 Hz

therefore I use a bandpass filter centered around 0.5 Hz

then there is a method to select the longest contiguous segment where we can identify this tone and not other frequencies

hope it helps !

NB : you will need the attached m function :

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

filename = 'Data.mat'; 
load(filename);
signal = data(:); % let's put the data in column oriented direction

Fs = 40; % sampling rate is 40 Hz here
dt = 1/Fs;

[samples,channels] = size(signal);
time = (0:samples-1)'*dt;

channels_selected = 1;

%% band pass  filter section %%%%%%
%bandpass filter
f_low = 0.25;
f_high = 1;
N_bpf = 2;           % filter order
[b,a] = butter(N_bpf,2/Fs*[f_low f_high]);
signal_filtered = filtfilt(b,a,signal);

       
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FFT parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
NFFT = 500;    % 2 s of data => 2 x Fs samples
OVERLAP = 0.95; % (75% ) overlap 
% spectrogram dB scale
spectrogram_dB_scale = 80;  % dB range scale (means , the lowest displayed level is XX dB below the max level)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display 1 : time domain plot
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

delta_signal = abs(signal - signal_filtered);
signal_rectified = abs(signal);
[delta_signal_env] = env_secant(time, delta_signal,30, 'top');
[signal_rectified_env] = env_secant(time, signal_rectified,30, 'top');

% logical statement
threshold = 0.04;
ind = (signal_rectified_env>threshold) & (delta_signal_env<threshold);
% define start / end points of contiguous segments
[begin,ends] = find_start_end_group(ind);
durations = ends-begin;

% select max duration segment (there will be only one selected)
[v,ind2] = max(durations);
ind = (begin(ind2):ends(ind2)) ;
time_selected = time(ind);
signal_selected = signal(ind);

figure(1),
subplot(2,1,1),plot(time,signal,'b',time,signal_filtered,'r',time_selected,signal_selected,'g');grid on
title(['Data Time plot  / Fs = ' num2str(Fs) ' Hz ']);
xlabel('Time (s)');ylabel('Amplitude');
legend('raw','filtered','selection')
subplot(2,1,2),plot(time,signal_rectified_env,'b',time,delta_signal_env,'r',time,threshold*ones(size(time)),'g--');grid on
title(['Data Time plot  / Fs = ' num2str(Fs) ' Hz ']);
xlabel('Time (s)');ylabel('Amplitude');
legend('signal rectified','delta signal','threshold')


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

[freq, spectrum_raw] = myfft_peak(signal,Fs,NFFT,OVERLAP);
[freq, spectrum_filtered] = myfft_peak(signal_filtered,Fs,NFFT,OVERLAP);
df = freq(2)-freq(1); % frequency resolution 

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display 3 : time / frequency analysis : spectrogram 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for ck = 1:channels
    [sg,fsg,tsg] = specgram(signal(:,ck),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

    % 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+ck);
    imagesc(tsg,fsg,sg_dBpeak);colormap('jet');
    axis('xy');colorbar('vert');grid on
    df = fsg(2)-fsg(1); % freq resolution 
    title([' Spectrogram / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(df,3) ' Hz / Channel : ' num2str(ck)]);
    xlabel('Time (s)');ylabel('Frequency (Hz)');

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [begin,ends] = find_start_end_group(ind)
% This locates the beginning /ending points of data groups
% Important : ind must be a LOGICAL array

    D = diff([0;ind(:);0]);
    begin = find(D == 1);
    ends = find(D == -1) - 1;
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 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

[CM Sahu]

Thank you for addressing my query.

Indeed, you're correct that the frequency of the additional sinusoidal wave is approximately 0.41Hz.

The approach is proving effective, although I'm still observing remnants of the sine wave in the filtered signal. Ideally, the designated region (highlighted in green) should register as zero following the elimination of the sine wave from the raw data. The range of superimposed sine wave is between the 613th and 1569th.

Please suggest any other method if you have.

Thank you very much

[Mathieu NOE]

hello

I tried to refine the code - now we can extract the sine signal in the 613th to 1569th sample window

and if you want to remove the sine wave from the total signal in that particular time window , this is also done in the code below :

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

filename = 'Data.mat'; 
load(filename);
signal = data(:); % let's put the data in column oriented direction

Fs = 40; % sampling rate is 40 Hz here
dt = 1/Fs;

[samples,channels] = size(signal);
time = (0:samples-1)'*dt;

channels_selected = 1;

ts = [613 1569]/Fs % this should be the time range where the sine signal must be detected and zeroed


%% band pass  filter section %%%%%%
%bandpass filter
f_low = 0.25;
f_high = 1;
N_bpf = 2;           % filter order
[b,a] = butter(N_bpf,2/Fs*[f_low f_high]);
signal_filtered = filtfilt(b,a,signal);

       
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FFT parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
NFFT = 500;    % 2 s of data => 2 x Fs samples
OVERLAP = 0.95; % (75% ) overlap 
% spectrogram dB scale
spectrogram_dB_scale = 80;  % dB range scale (means , the lowest displayed level is XX dB below the max level)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display 1 : time domain plot
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

signal_rectified = abs(signal);
[signal_rectified_env] = env_secant(time, signal_rectified,100, 'top');

signal_filtered_rectified = abs(signal_filtered);
[signal_filtered_rectified_env] = env_secant(time, signal_filtered_rectified,100, 'top');

% logical statement
threshold = 0.03;
ind = (signal_filtered_rectified_env>threshold);
% define start / end points of contiguous segments
[begin,ends] = find_start_end_group(ind);
durations = ends-begin;

% select max duration segment (there will be only one selected)
[v,ind2] = max(durations);
ind = (begin(ind2):ends(ind2)) ;
time_selected = time(ind);
signal_selected = signal(ind);

% refinement of the detection of the beginning and end of sine signal 
tmp = abs(diff(abs(signal_selected)));
[p,locs] = findpeaks(tmp,'MinPeakHeight',max(tmp)/10);
l1 = locs(1)-1;
l2 = locs(end)+1;
ind = ind(l1:l2);
time_selected = time(ind);
signal_selected = signal(ind);

figure(1),
subplot(2,1,1),plot(time,signal,'b',time,signal_filtered,'r',time_selected,signal_selected,'g');grid on
title(['Data Time plot  / Fs = ' num2str(Fs) ' Hz ']);
xlabel('Time (s)');ylabel('Amplitude');
legend('raw','filtered','selection')
subplot(2,1,2),plot(time,signal_rectified_env,'b',time,signal_filtered_rectified_env,'r',time,threshold*ones(size(time)),'g--');grid on
title(['Data Time plot  / Fs = ' num2str(Fs) ' Hz ']);
xlabel('Time (s)');ylabel('Amplitude');
legend('signal rectified','signal filtered rectified','threshold')

% remove the selected portion (zero it)
signal_corrected = signal;
signal_corrected(ind) = 0;

figure(2),
plot(time,signal,'b',time,signal_corrected,'r');grid on
title(['Data Time plot  / Fs = ' num2str(Fs) ' Hz ']);
xlabel('Time (s)');ylabel('Amplitude');
legend('raw','corrected')

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

[freq, spectrum_raw] = myfft_peak(signal,Fs,NFFT,OVERLAP);
[freq, spectrum_filtered] = myfft_peak(signal_filtered,Fs,NFFT,OVERLAP);
df = freq(2)-freq(1); % frequency resolution 

figure(3),
semilogx(freq,20*log10(spectrum_raw),'b',freq,20*log10(spectrum_filtered),'r');grid on
title(['Raw data Averaged FFT Spectrum  / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(df,3) ' Hz ']);
xlabel('Frequency (Hz)');ylabel('Amplitude (dB)');
legend('raw','filtered')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [begin,ends] = find_start_end_group(ind)
% This locates the beginning /ending points of data groups
% Important : ind must be a LOGICAL array

    D = diff([0;ind(:);0]);
    begin = find(D == 1);
    ends = find(D == -1) - 1;
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 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

[CM Sahu]

Thank you sincerely for your dedication in tackling this issue.

However, it appears that the selected segment (highlighted in green) is inadvertently removing everything, including the underlying signal upon which this sinusoidal wave is overlaid.

Your efforts are greatly appreciated, but could you please endeavor to isolate only the sine wave component from the original signal? If I were to express this mathematically, it would be represented as follows:

y(t) = A1*sin(2*pi*f*t) + x(t),

Here, x(t) denotes the original signal, and A1*sin(2*pi*f*t) represents the superimposed wave between the 613th and 1569th data points. Our objective is to solely remove this portion A1*sin(2*pi*f*t) from the modified signal y(t). Thank you for your attention to this matter.

[Mathieu NOE]

Ok, so this is important info , the sine wave component is supposed to be a constant sinus fonction from t = 613th to 1569th sample ?

see the amplitude is dropping after t = 20s (800th sample) , so I were to fit a fixed sinewave model , I would then see a sinoisidal error showing up at the same frequency and you would probably say that that shoudn't be there

I think my approach of simply bandpass filtering is not perfectl because it will extract more than just one sinewave and of course if you remove more than the sinewave it's not going to give you the expected result - just to show what we would get , it's here :

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

filename = 'Data.mat'; 
load(filename);
signal = data(:); % let's put the data in column oriented direction

Fs = 40; % sampling rate is 40 Hz here
dt = 1/Fs;

[samples,channels] = size(signal);
time = (0:samples-1)'*dt;

channels_selected = 1;

ts = [613 1569]/Fs % this should be the time range where the sine signal must be detected and zeroed


%% band pass  filter section %%%%%%
%bandpass filter
f_low = 0.25;
f_high = 1;
N_bpf = 2;           % filter order
[b,a] = butter(N_bpf,2/Fs*[f_low f_high]);
signal_filtered = filtfilt(b,a,signal);

       
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FFT parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
NFFT = 500;    % 2 s of data => 2 x Fs samples
OVERLAP = 0.95; % (75% ) overlap 
% spectrogram dB scale
spectrogram_dB_scale = 80;  % dB range scale (means , the lowest displayed level is XX dB below the max level)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display 1 : time domain plot
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

signal_rectified = abs(signal);
[signal_rectified_env] = env_secant(time, signal_rectified,100, 'top');

signal_filtered_rectified = abs(signal_filtered);
[signal_filtered_rectified_env] = env_secant(time, signal_filtered_rectified,100, 'top');

% logical statement
threshold = 0.03;
ind = (signal_filtered_rectified_env>threshold);
% define start / end points of contiguous segments
[begin,ends] = find_start_end_group(ind);
durations = ends-begin;

% select max duration segment (there will be only one selected)
[v,ind2] = max(durations);
ind = (begin(ind2):ends(ind2)) ;
time_selected = time(ind);
signal_selected = signal(ind);

% refinement of the detection of the beginning and end of sine signal 
tmp = abs(diff(abs(signal_selected)));
[p,locs] = findpeaks(tmp,'MinPeakHeight',max(tmp)/10);
l1 = locs(1)-1;
l2 = locs(end)+1;
ind = ind(l1:l2);
time_selected = time(ind);
signal_selected = signal(ind);

figure(1),
subplot(2,1,1),plot(time,signal,'b',time,signal_filtered,'r',time_selected,signal_selected,'g');grid on
title(['Data Time plot  / Fs = ' num2str(Fs) ' Hz ']);
xlabel('Time (s)');ylabel('Amplitude');
legend('raw','filtered','selection')
subplot(2,1,2),plot(time,signal_rectified_env,'b',time,signal_filtered_rectified_env,'r',time,threshold*ones(size(time)),'g--');grid on
title(['Data Time plot  / Fs = ' num2str(Fs) ' Hz ']);
xlabel('Time (s)');ylabel('Amplitude');
legend('signal rectified','signal filtered rectified','threshold')

% remove the sine wave in the selected portion 
signal_corrected = signal;
signal_corrected(ind) = signal(ind) - signal_filtered(ind);

figure(2),
plot(time,signal,'b',time,signal_corrected,'r');grid on
title(['Data Time plot  / Fs = ' num2str(Fs) ' Hz ']);
xlabel('Time (s)');ylabel('Amplitude');
legend('raw','corrected')

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

[freq, spectrum_raw] = myfft_peak(signal,Fs,NFFT,OVERLAP);
[freq, spectrum_filtered] = myfft_peak(signal_filtered,Fs,NFFT,OVERLAP);
df = freq(2)-freq(1); % frequency resolution 

figure(3),
semilogx(freq,20*log10(spectrum_raw),'b',freq,20*log10(spectrum_filtered),'r');grid on
title(['Raw data Averaged FFT Spectrum  / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(df,3) ' Hz ']);
xlabel('Frequency (Hz)');ylabel('Amplitude (dB)');
legend('raw','filtered')

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display 3 : time / frequency analysis : spectrogram 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for ck = 1:channels
    [sg,fsg,tsg] = specgram(signal(:,ck),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

    % 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+ck);
    imagesc(tsg,fsg,sg_dBpeak);colormap('jet');
    axis('xy');colorbar('vert');grid on
    df = fsg(2)-fsg(1); % freq resolution 
    title([' Spectrogram / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(df,3) ' Hz / Channel : ' num2str(ck)]);
    xlabel('Time (s)');ylabel('Frequency (Hz)');

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [begin,ends] = find_start_end_group(ind)
% This locates the beginning /ending points of data groups
% Important : ind must be a LOGICAL array

    D = diff([0;ind(:);0]);
    begin = find(D == 1);
    ends = find(D == -1) - 1;
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 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

[CM Sahu]

Thank you sincerely for your hard work.

I believe it's impossible to isolate this sine wave from the original signal. The challenge lies in the sine wave's fluctuating amplitude and frequency. Despite attempting to sweep the frequency of the sine wave over time, unfortunately, I wasn't successful.

Regarding this enclosed figure:

- The red line represents the original signal graph.

- The blue line indicates the superimposed sinusoidal wave on the original wave.

- The green line illustrates the graph generated by your suggested method.

Thank you for your effort!

[Mathieu NOE]

tx for the nice comment

I see between the red and blue that the sinewave is not stationnary

so the model you suggested y(t) = A1*sin(2*pi*f*t) + x(t), is not viable, unless you accept larger error because the "sine" is (again) unstationnary

so it's quite challenging to decide what belongs to the sine wave and what is the original signal

[CM Sahu]

Ok, Thank you for your effort

[Mathieu NOE]

as always, my pleasure !