hello
I notice that the cutoff freq is very low compared to fs , so maybe lowpass has a numerical problem when the normalized cut off freq (f_cutoff / fs) is close to zero
I plotted some sections of the (time domain) signals and it seems indeed there is no filtering effect to be seen (with high fs)
more tests would be needed to double check what's going on with lowpass function
I tried also to reduce fs by a factor of 10 and then it seemed to work better, but that is maybe not what you want (to reduce fs)
so instead I would suggest to come back to the "good old" way to design your filter combined with filtfilt , and then it seems to work just fine even with the high fs value;
I also put the x scale in log for the fft spectrum plots as it shows more clearly that the spectrum below f_cutoff remains the same , and the filtering above is now effective
hope it helps
%Parameter Setting
pixel_image = 256; %Input the number of pixels in the image obtained by raster scanning (input 2^n)
dr = 1/(2*sqrt(3)); %Enter dither circle radius [grid].
a_fast_grid = 10; %fast axis scanning range [grid]
a_slow_grid = 10; %Slow axis scanning range [grid]
fm=5000; %Dither circle modulation frequency [Hz]
fs= fm*240 ; %Sampling frequency [Hz]
f_fast = 10.2; %Input scanning frequency [Hz] (1 line scanning count in 1[s])
start_point_x = 0; %Input x-coordinate of scanning start point (input 1 if you want to move by 1[grid])
start_point_y = 0; %Input y-coordinate of scanning start point (input 1 if you want to move by 1[grid])
%Parameter setting for fast-axis triangular wave
amplitude_fast = a_fast_grid/2; %fast axis amplitude
%Parameter setting for slow-axis triangular wave
amplitude_slow = a_slow_grid/2; %slow axis amplitude
f_slow = (f_fast)/(2*pixel_image); %Slow axis triangular wave frequency
% Generation of time vectors
total_time=256/f_fast; %Total Scan Time
t = linspace(0, total_time, fs * total_time);
x_raster = start_point_x + amplitude_fast*(2/pi)*acos(cos(2*pi*f_fast*t));
y_raster = start_point_y + amplitude_slow*(2/pi)*acos(cos(2*pi*f_slow*t));
x_dither = dr*cos(2*pi*fm*t);
y_dither = dr*sin(2*pi*fm*t);
x = x_raster + x_dither;
y = y_raster + y_dither;
z1 = cos(2*pi*((x-y)/(sqrt(3))));
z2 = cos(2*pi*(2*y/(sqrt(3))));
z3 = cos(2*pi*((x+y)/(sqrt(3))));
z = (z1 + z2 + z3);
%Reference Signal Creation
w = 2*pi*fm ;
phi = 0 ;
reference_signal = sin(w*t+phi) ;
%Mixing signal creation
mixising_signal = z.* reference_signal ;
%Applying a low-pass filter
f_cutoff = 500 ; %Cutoff frequency [Hz] setting
% lowpass_signal = lowpass(mixising_signal,f_cutoff,fs) ;
%% alternative solution :
[zer,pol,k] = butter(6,f_cutoff/(fs/2));
[SOS,G] = zp2sos(zer,pol,k);
lowpass_signal = filtfilt(SOS, G, mixising_signal);
% NB : It is best not to use the transfer function implementation of discrete (digital) filters,
% because those are frequently unstable.
% Use zero-pole-gain and second-order-section implementation instead.
%The result of FFT.-------------------------------------------
n = length(t); % Number of samples
f = (0:n-1)*(fs/n); % Frequency axis [Hz]
half_n = floor(n/2); % Half of data
fft_z = abs(fft(z) / n);
fft_mixising_signal = abs(fft(mixising_signal) / n);
fft_lowpass_signal = abs(fft(lowpass_signal) / n);
figure;
tiledlayout(3,1)
nexttile
semilogx(f(1:half_n), fft_z(1:half_n));
title('FFT of tunnel current'); xlabel('Frequency [Hz]'); ylabel('Amplitude');
xlim([0 20000]),ylim([0 0.3]);
nexttile
semilogx(f(1:half_n), fft_mixising_signal(1:half_n));
title('FFT of mixing signal'); xlabel('Frequency [Hz]'); ylabel('Amplitude');
xlim([0 20000]),ylim([0 0.3]);
nexttile
semilogx(f(1:half_n), fft_lowpass_signal(1:half_n));
title('FFT after applying low-pass'); xlabel('Frequency [Hz]'); ylabel('Amplitude');
xlim([0 20000]),ylim([0 0.3]);
