%% 1-D signal

AC=255/2;

DC=255/2;

f0=4; %frequency

ph=pi/6; %phase offset

x=linspace(0,1,1024);

ysq=DC+AC*square(2*pi*f0*x+ph);

figure; plot(x,ysq); title('square wave');

% Fourier analysis

sq_fft=fft(ysq);

mg_sq=abs(sq_fft);

mg_sq=mg_sq/length(mg_sq);

P=1+ceil((length(mg_sq)-1)/2); %Position of the maximum frequency

mg_sq(2:P)=2*mg_sq(2:P);

mg_sq=mg_sq(1:P); %Keep only positive frequencies

figure; stem(0:length(mg_sq)-1,mg_sq); title('magnitude spectrum')

%Retrieval of individient components by keep only the positive frequency

%f=n*f0, n is odd;

%amp=(2*A)/(n*pi), A=255, n=odd

fft_f4=zeros(size(sq_fft)); %freq = 4;

fft_f4(5)=sq_fft(5);

fft_f12=zeros(size(sq_fft)); %freq = 4;

fft_f12(13)=sq_fft(13);

figure; plot(x,2*real(ifft(fft_f4)),'r-'); hold on;

plot(x,2*real(ifft(fft_f12)),'b-');

title('Retrieved components');

%Alternatively: keep both negative and positive frequncies

fft_f4=zeros(size(sq_fft)); %freq = 4;

fft_f4(5)=sq_fft(5); fft_f4(end-3)=sq_fft(end-3);

fft_f12=zeros(size(sq_fft)); %freq = 4;

fft_f12(13)=sq_fft(13); fft_f12(end-11)=sq_fft(end-11);

figure; plot(x,real(ifft(fft_f4)),'r-'); hold on;

plot(x,real(ifft(fft_f12)),'b-');

title('Retrieved components');

%or shifting and anti-shifting the frequency origin

sq_sfft=fftshift(sq_fft); % The central position: P, 0 frequency

fft_f4=zeros(size(sq_fft));

fft_f4(P-4)=sq_sfft(P-4); fft_f4(P+4)=sq_sfft(P+4); fft_f4=fftshift(fft_f4);

fft_f12=zeros(size(sq_fft));

fft_f12(P-12)=sq_sfft(P-12); fft_f12(P+12)=sq_sfft(P+12); fft_f12=fftshift(fft_f12);

figure; plot(x,real(ifft(fft_f4)),'r-'); hold on;

plot(x,real(ifft(fft_f12)),'b-');

title('Retrieved components');