There are a couple issues with your code. First, let's talk about why you can do what you're trying to do.
What you're trying to take advantage of is the duality between a multiplication in the time domain and convolution in the frequency domain. The fftshift() function simply swaps the two halves of a vector (which is assumed to be the result of a call to fft()). This operation is equivalent to a circular rotation of half of the vector length in the frequency domain.
One can accomplish the same operation by multiplying by a complex exponential function in the time domain before the FFT, as you noted. To get an $M$-sample shift, you need to multiply by a complex exponential that looks like:
$$
s[n] = e^{\frac{j2\pi nM}{N} }
$$
For the case where $M = \frac{N}{2}$, this simplifies to:
$$
s[n] = e^{j\pi n} = [1, -1, 1, -1, \ldots]
$$
This illustrates problem #1 with your script: it multiplies by the inverse of the above: $-s[n]$. You need to invert the sense of the if statement in the for loop to ensure that the vector s has the appropriate value.
Problem #2 stems from another simple MATLAB mistake: all of the vectors in your simulation are length 50001. So, the time-domain rotation that you're doing by multiplying by the $s[n]$ illustrated above actually results in an effective frequency domain shift of 25000.5 bins (instead of the 25000-bin shift implemented by fftshift()). The mismatch in results is due to this extra half-bin shift.
You can fix this problem in two ways:
Change the time-domain multiplication to multiply by an actual complex exponential vector at the appropriate frequency (i.e. set $M=25000$ and $N=50001$ in the first equation). The resulting signal at the input to the FFT will be complex.
Change your simulation to use even-length vectors (i.e. chop off the last element of the t vector that you initialize at the beginning).
I tweaked your simulation to fix the generation of s and truncate the vector length to 50000. It is as follows:
%1khz
f=1000;
t=0:.0001:5;
t=t(1:length(t)-1);
x=sin(2*pi*f*t);
%Wave
subplot(4,1,1);
plot(x)
%FFT
subplot(4,1,2)
ffted=fft(x);
plot(real(ffted));
%FFT Shifted
subplot(4,1,3)
shifted=fftshift(ffted);
plot(real(shifted));
%FFT Shifted by preshifting
subplot(4,1,4)
s=zeros(1,length(x));%In languages like C++ this is computationally inexpesnive
for i=1:length(x)
if (~bitget(i,1))%odd evens, works for real signal
s(i)=-1;
else
s(i)=1;
end
end
chirpshifted=fft(s.*x);
plot(real(chirpshifted));
I then compare the results of the two methods:
>> diff = abs(chirpshifted-shifted);
>> sum(diff)
ans =
3.3406e-21
As you can see, the results agree as perfectly as you would expect from floating-point math.
