The tool of choice here is the z-transform, which is applied to your equation. It transforms sequences of complex numbers to functions in the complex plane. Its most useful properties are that it takes delays to powers of the argument z, and for
z = exp(i omega) we get the fourier transform of the transformed sequence. Look up the details, I'll just show you how it works here.
y[n] = x[n] - x[n-486] ---> Y(z) = X(z) - X(z)*z^(-486)
The equation in z-domain factors on the right hand side to give
Y(z) = X(z) ( 1 - z^(-486) ).
Since Y(z) is the output and X(z) is the input, for linear systems we define the transfer function to be
Y(z)/X(z) = H(z) and the above equation becomes
H(z) = Y(z)/X(z) = 1 - z^(-486)
Now if we substitute z = exp(i omega) we get the frequency response of your system, which is
H_f(omega) = 1 - exp(i omega)^(-486)
and with the properties of exp() it simplifies to
H_f(omega) = 1 - exp(- 486 i omega)
If you need actual frequency f instead angular frequency omega, just substitute
omega = 2 pi f/fs, with sampling rate fs.
Finally, you'll probably want to know just the magnitude of the response, so look at abs(H_f).
In this specific case it's very simple to calculate because
abs(H_f) = sqrt( H_f * conj(H_f) ) and
conj(H_f) = 1 - exp(+ 486 i omega).
As requested, some code. Here's what you can do in matlab/octave to get the magnitude frequency response with a logarithmic frequency axis and magnitude axis in dB.
% Create a logarithmically spaced frequency axis with range 20Hz to 20kHz
fAxis = logspace( log10( 20 ) , log10( 20000) , 10000);
% calculate the complex frequency response at these frequencies
% Set SampleRate before calling this, like SampleRate = 44100;
Hf = 1 - exp(1i*486*2*pi*fAxis/SampleRate);
% plot the magnitude spectrum