matlab fft phase response square wave

Question

I am new to phase analysis, recently I have been trying compare the phase between the input signal and output signal of a system. From the code below, I simulated a square wave and plotted the phase, from basic theory, I thought that since this is a simulated signal, I should get near zero values other than the values at each input frequency. So my questions are (based on the following code):

1) Why does Matlab give me pi/2 as phase response of actual input signals, I don't understand the calculation behind it (does it have to do with my signal being sine and not cosine?)

2) If I don't filter out low (near zero) values, the phase response is a sloping line that goes toward very large numbers, and seems to be proportional to the number of samples of FFT I have, is this due to how unwrap works?

3) If I am trying to compare an input signal vs output signal, both of which for the sake of argument looks similar to the square wave I simulated, how do I determine if the system is linear phased? (from theory, the phase delay should be the same therefore, after finding the phase response, should I check if the difference of phase for each input signal between output and input is constant?

4) Also, group delay is considered as the derivative of phase response, and phase delay is the normalization of phase with respect to its corresponding frequency, a perfect square wave should have constant values at each harmonic for group delay and phase delay shouldn't it?, if so can I somehow see this result in my simulation?

close all;

clear all;

fs = 16000;

end_time = 6-1/fs;

x = 0:1/fs:end_time;

x = x';

n=2000;

freq = 32;

wave = sin(2*pi* freq*x)/(4*pi);

for i = 3:2:n

    wave = wave + sin(i*2*pi*freq*x)/(4*i*pi);

end

wave = awgn(wave,55,'measured');

wave (1:10000) = 0;

NFFT = length(wave);

fft_result = fft(wave, NFFT);


for i = 1:NFFT

    if abs(fft_result(i)) < 10^-4

        fft_result(i) = 0;

    end

end

fft_mag = db(abs(fft_result));

fft_phase = unwrap(angle(fft_result));

fstep = fs/NFFT;

f_range = 0:fs/NFFT:fs-1/NFFT;


figure

plot(f_range, fft_phase);

Thank you for being patient enough to read all of this and giving me answers. =)

Answer 1

Re question (1):

Most FFTs, or the way atan2() is commonly used with complex FFT results to compute a phase, produce a zero phase (strictly real) output for strictly symmetric input. This means you want cosine waves, periodic in aperture, not sine waves starting at zero at the window edge, for FFT results with zero phase.