My butterworth filter isn't preventing aliasing

Hi all,

I recently read this guide to aliasing here:

http://www.mstarlabs.com/dsp/antialiasing/antial.html

And I wanted to make a demonstration in MATLAB. What I want to demonstrate is this: high frequency signals that are above your sample rate will alias onto lower frequencies. This aliasing can be prevented by low-pass filtering the signal before it's sampled.

However, a Fourier transform of my filtered signal still shows an alias. Could this be just because my butterworth filter is a digital filter, rather than analog?

Below is my code, there's signals at 50 and 60 Hz, and then a high frequency signal at 11,000 Hz. The sample rate is 10,000 Hz, and the cutoff frequency is 2,500 Hz. I have added noise, as well, to demonstrate that the filter is working at the right frequency. The 11,000 Hz signal is being aliased onto 1,000 Hz, when it should be eliminated.

Thanks in advance!

clear;
close all;
clc;

%__________________________________________________________________________
% Inputs
Fs =10000;                    % Sampling frequency
Max_time = 10;                % Signal duration in sec
w_a = 30;                     % Signal A frequency
w_b = 75;                     % Signal B frequency
w_c = 11000;                  % Signal C frequency
noise_gain = 10;              % Gain applied to noise
cutoff_freq = 2500;

%__________________________________________________________________________
% Calculations
T = 1/Fs;                     % Sample time
L = Fs*Max_time;              % Number of samples
t = (0:L-1)*T;                % Time vector

cutoff_rad_persec = cutoff_freq*2;
cutoff_rad_sample = cutoff_rad_persec/Fs;

% Sum of 3 signals
x = sin(2*pi*w_a*t) + sin(2*pi*w_b*t) + sin(2*pi*w_c*t); 
y = x + noise_gain*randn(size(t));     % Sinusoids plus noise

% Filter the output
[b,a] = butter(5,cutoff_rad_sample);
y_filtered = filter(b,a,y);

% Fourier transform of signal
NFFT = 2^nextpow2(L);         % Next power of 2 from length of y
Y = fft(y,NFFT)/L;
Y_filtered = fft(y_filtered,NFFT)/L;

% Make a frequency axis ranging from 0 to the Nyquist frequency, 500 Hz (Fs/2)
f = Fs/2*linspace(0,1,NFFT/2+1);

%__________________________________________________________________________
% Plot data
FigHandle = figure;
set(FigHandle, 'Position', [0, 40, 2000, 1070]);
set(0,'DefaultAxesFontSize', 22)
set(0,'defaultlinelinewidth',1)

% Plot in time domain
subplot(2,1,1);
plot(t,y)
hold on
plot(t,x,...
    'LineWidth',1.5)

title([...
    'Sum of 3 Signals (',...
    num2str(w_a),...
    ' Hz, ',...
    num2str(w_b),...
    ' Hz, and ',...
    num2str(w_c),...
    ' Hz) Corrupted with Zero-Mean Random Noise'])
xlabel('Time (s)')
xlim([0 .16])
ylim([-4 4])
box off

% Plot in frequency domain, single-sided amplitude spectrum.
subplot(2,1,2);

% Plot of unfiltered signal
%plot(f,2*abs(Y(1:NFFT/2+1)))
%hold on

% Plot of filtered signal
plot(f,2*abs(Y_filtered(1:NFFT/2+1)))
hold on

% Plot a line at the cutoff frequency
line_x = [cutoff_freq cutoff_freq];
line_y = [0 abs(max(2*Y))*1.1];
plot(line_x, line_y,...
    'g-',...
    'LineWidth',1.2);

title('Frequency Analysis, Single-Sided Amplitude Spectrum of y(t)')
xlabel('Frequency (Hz)')
ylabel('|Y(f)|')
ylim([0 1])
box off