As i understand you are failing to understand the output phase associated with sinusoid of frequency 4000 Hz.
Let us look at the following equation
s_1 = cos(2*pi*1500*t1 - deg2rad(20)) - sin(2*pi*4000*t1 - deg2rad(55));
the above can be rewritten as
s_1 = cos(2*pi*1500*t1 - deg2rad(20)) + cos(2*pi*4000*t1 + deg2rad(90) - deg2rad(55));
we have used the property cos(90 + x) = -sin(x) so it simplifies to the following
s_1 = cos(2*pi*1500*t1 - deg2rad(20)) + cos(2*pi*4000*t1 + deg2rad(35));
the number 35 is a result of simple angle addition of 90 and -55.
So in the phase spectra you will see the phase 35 degrees associated with sinusoid of 4000 Hz.
The uploaded code does not run properly unless few modifications are made, here is a revised code snippet
Fs = 400e3; % Sampling rate
T = 1/Fs ; % Sampling period
N = 40000; % Signal length
t1 = (0:N-1)*T % Time vector
s_1 = cos(2*pi*1500*t1 - deg2rad(20)) - sin(2*pi*4000*t1 - deg2rad(50));
%Getting FFT
y1 = fft(s_1);
z1 = fftshift(y1);
l = length(y1);
figure ;
f1 = (-l/2:l/2-1)/l*Fs;
stem(f1,abs(z1))
xlabel('Frequency (Hz)')
ylabel('|y1|')
grid
%Getting phase
tol1 = 1e-6;
z1(abs(z1) < tol1) = 0;
theta1 = rad2deg(angle(z1));
figure;
stem(f1,theta1);
xlabel('Frequency (Hz)')
ylabel('Phase / \pi')
grid
I hope this answers your queries, please accept the query as answered if satisfied !
