You seem to be confused with the time and frequency domain. Your code for an LP filter is the inverse Fourier transform of an Ideal low-pass filter. So by getting a frequency response for the function you wrote you'll get an ideal LP response you want.
However, there are two mistakes:
1) Your frequency for the ideal LP filter needs to be $[-\frac{m-1}{2} : \frac{m-1}{2}]$
2) You need to change the value for the 0-frequency only if M is an odd number.
Now, you can use the function freqz
to get your freqency response, so your final code would be like this:
%[y,Fs]=wavread('sound.wav');
%m = length(y); %30586
m = 30586;
Fc = 500; %% frequency cut
n = -(m-1)/2 : (m-1)/2; %% kernel length / samples
%Fs = Fs; %22000 %% sample frequency / (nyquist - %2*m)
Fs = 22000;
Tc = 2*pi*(Fc/Fs); %% theta c
Hm = sin(Tc*n) ./ (pi*n);
if (rem(m, 2) == 1)
index = (m+1)/2;
Hm(index) = Tc/pi;
end
[H w] = freqz(Hm, 1, 1024); % 1024 - or how many points you want for your Fourier Transform
% Now your H is your frequency response, and w is digital frequency in [0,pi].
% You can now plot the amplitude response:
figure,plot(w, abs(H));
And you have your Low-Pass filter. However, I'd like to suggest using the many ways to create filters in Matlab, they were made for a reason, and allow for much easier filter design :)
One more remark, I have to suggest that you know the difference between analog and digital frequency, as they have different notation. Analog frequencies you write in uppercase, e.g. $F, \Omega$, and digital frequencies in lower-case, $f, \omega$. Here's a version of the code which is a little prettier:
M = 30586;
Fc = 500;
Fs = 22000;
wc = 2*pi*Fc/Fs;
n = -(M-1)/2:(M-1)/2;
h = (sin(n*wc))./(n*pi);
if ( M-2*fix(M/2) ) > 0
index = (M+1)/2;
h(index) = wc/pi;
display(index);
end
[H w] = freqz(h, 1, 1024);
plot(w, abs(H));
Cheers :)
EDIT:
When you use freqz, the digital frequency w can be somewhat hard to control, because it's not created in the same way, or more clearly, the points are not selected like a normal frequency we create with, for example, F = 0 : Fs/(M-1) : Fs
, even though it should be the "same" frequency.
The result is that it's not really possible to get the output filtered signal with Y = H .* X
, and y = ifft(Y)
.
However, I would suggest using the function filter(b, a, x)
for filtering your signal, which works good with freqz
. You can proceed to do just that:
% ------
% Create a random signal x
n = 0:M-1;
F = -Fs/2:Fs/(M-1):Fs/2;
Fx = [2000 7000 10000];
x = cos(2*pi*Fx(1)/Fs*n) + cos(2*pi*Fx(2)/Fs*n) + cos(2*pi*Fx(3)/Fs*n);
X = abs(fftshift(fft(x, Nfft))); % fftshift because we need domain [-Fs/2, Fs/2]
% ---------
% -------
% Now filter the signal
y = filter(H, 1, x);
Y = abs(fftshift(fft(y)));
% -------
That's it :)