How do I remove frequency component on my spectrum after generating interference signal?

Hi, I simulated a interference signal which is based on Michelson's interferometer.

However, in real case, there is a frequency component on my spectrum. (Figure 2)

This is a perfect spectrum.

This has a frequency component on Figure 1.

After doing FFT, interference signal was generated.

The signal at about 57 kHz is the main signal. Besides this main signal, the other two peaks are because of frequency component on the original spectrum.

Therefore, how can I suppress this two signal which is caused by frequency component after generating interference signal?

clear all
clf
close all
clc


%% commercial light source
filename = ('APD_commercial_spectrum.xlsx');
data = xlsread(filename);
points = data(11:1120,1);
load('wavelength.mat');
load('depth.mat');
signal = data(:,2);


figure
FFT = FFT_hann(signal,-100,-20,'PSF of commercial spectrum');

%%Time specifications:
Fs = 200e3;                   % samples per second
dt = 1/Fs;                   % seconds per sample
t = linspace(dt,6e-3,1200)';     % seconds


%% first one: shape of commercial light source
num = 1200;
lambda_s = 770*10^-9; %start wavelength
lambda_f = 830*10^-9; %end wavelength
k_s =2*pi/lambda_f;%start k-space
k_f =2*pi/lambda_s; %end k-space
k_c = (k_s+k_f)/2; %center k-space
k = linspace(k_s,k_f,num)';


mask = signal > 1000;
maskedK = k(mask);
maskedSignal = signal(mask);
coefficients = polyfit(maskedK, maskedSignal, 2);
% Get smoothed signal.
smoothedSignal = polyval(coefficients,k(11:1120));
% smoothedSignal = rescale(smoothedSignal);
empty = zeros(10,1);
empty2 = zeros(80,1);
S = [empty;smoothedSignal;empty2];
profile = rescale(S);
commercial = profile;
k = linspace(k_s,k_f,length(commercial))';

%% plot spectrum
figure
plot(t,commercial,'linewidth',1.5);
xlabel('Time (s)');
ylabel('Amplitude');
title('Simulate a commercial spectrum');
set(gca,'fontsize',14,'linewidth',1.5);


figure
plot(wavelength,Rs,'linewidth',1.5);
title('Reflectivity of sample arm');
xlabel('wavelength (nm)');
ylabel('Reflectivity');
zoom xon;
set(gca,'fontsize',14,'linewidth',1.5);

%%
%frequency on envelope
noise = 0.05*rand(size(profile));
freqs = 0.03.*cos(2.*k.*369*1e-6);

S_f = commercial+freqs + noise;
figure
plot(t,S_f,'linewidth',1.5);
xlabel('Time (s)');
ylabel('Amplitude');
title('spectrum with frequency components and noise');
set(gca,'fontsize',14,'linewidth',1.5);

figure
FFT = FFT_hann(S_f,-100,-20,'Spectrum with frequency components and noise');

%%

interference = IN(S_f,2000);
% interference = rescale(interference);
figure
FFT = FFT_hann(interference,-100,-50,'Interference signal');
hold on


function interference = IN(S_f,depth);

load('R_770_825.mat');
load('coupler_770_825.mat');

Rm_ref = R./100;
Rm_sam = R./100;

%%Time specifications:
Fs = 200e3;                   % samples per second
dt = 1/Fs;                   % seconds per sample
StopTime = 0.25;             % seconds
t = linspace(dt,6e-3,1200)';     % seconds

ratio1 = y1./100; %coupling ratio
ratio2 = y2./100;

num = 1200;
lambda_s = 770*10^-9; %start wavelength
lambda_f = 825*10^-9; %end wavelength
k_s =2*pi/lambda_f;%start k-space
k_f =2*pi/lambda_s; %end k-space
k_c = (k_s+k_f)/2; %center k-space
k = linspace(k_s,k_f,num)'; 

I_ref = ratio1.*ratio2.*Rm_ref;
I_sam = ratio1.*ratio2.*Rm_sam;

interference = S_f.*(2.*sqrt(I_sam.*I_ref).*(cos(2.*k.*depth*1e-6)));

end


function PSF = FFT_hann(signal,yscale,threshold,Titlename);
load('depth.mat');
load('wavelength.mat');

len = length(signal);
%hanning window length
win = hanning(len);
%multiply a hanning window for great resolution
signal_win = signal.*win;
%doing FFT and normalized
PSF = abs(fftshift(fft(signal_win,2048)));
PSF = PSF./max(PSF);

wav = linspace(770,830,len);

% figure
subplot(2,1,1)
plot(wav,signal,'linewidth',1);
title(sprintf('%s',Titlename));
% title('Original spectrum');
xlabel('wavelength (nm)');
ylabel('Amplitude');
set(gca,'fontsize',15,'linewidth',1);
grid on
hold on

logPSF = 10*log(PSF);


subplot(2,1,2)
plot(depth,logPSF,'linewidth',1);
% hold on
title('A-scan');
xlabel('Depth (µm)');
ylabel('Amplitude');
xlim([0 3500]);
ylim([yscale 0]);
set(gca,'fontsize',15,'linewidth',1);
grid on

%find peaks
[psor,lsor] = findpeaks(logPSF,'MinPeakHeight',threshold,'MinPeakDistance',40);
%round to two decimal places
depth=roundn(depth,0);
psor=roundn(psor,-2);

%if X=numel(lsor)/2 is odd number, rounds each element of X to the nearest integer greater than or equal to that element
 for t = ceil(numel(lsor)/2)+1:numel(lsor)
     %mark coordination 
%      text(depth(lsor(t)),psor(t)+5,['(',num2str(depth(lsor(t))),',',num2str(psor(t)),')'],'fontsize',12);
     text(depth(lsor(t))+15,psor(t),[num2str(psor(t)),' dB'],'fontsize',12);
 end

end

Answer 1

[Mathieu NOE]

hello

I have already answered this question here :

How to remove relative frequency information in time domain? - (mathworks.com)

all the best