# matlab How to design lpf\bpf\hpf without builtin functions

I'm trying to create a simple low pass filter but i'm not sure what is wrong.

my code:

[y,Fs]=wavread('sound.wav');
m = length(y); %30586
Fc = 500;               %% frequency cut
n = [0:m-1];            %% kernel length / samples
Fs = Fs;    %22000      %% sample frequency / (nyquist - %2*m)
Tc = 2*pi*(Fc/Fs);      %% theta c
Hm = sin(Tc*n) ./ (pi*n);
Hm(1) = Tc/pi;
figure,plot(Hm);


1.

LPF: How to get to the 0-500 range?

We used this code in the class but with Fc=45,m=100,Fs=200 and it looks ok. but when I try this on this wav file it's not working as expected. I thought it will pass from 0-500 but instead it returns some crowled plot.

only when I set Fc to be a low number I can see the sin wave. so i've try to manipulate by setting Fc=real(ceil(m/500)). it better than nothing but it's not it.

2.

hpf = 1-lpf, right? well is seems not to be. my guess is that my lpf is wrong

3.

I will deal the ripple after i'll manage to do this

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 :)

• thank you. the reason I did from 0-m is because I thought I need to shift the sin to zero. – jico Jun 2 '13 at 8:48
• maybe you too fast for me. I'm not sure how to apply the filter. I've set [H w] = freqz(h, 1, M); and then filtered like that: filtered=H .* fft(y); multiplication in the frequency domain. is that right? – jico Jun 2 '13 at 8:56
• it seems like i'm cutting Fc*2. – jico Jun 2 '13 at 9:06
• The problem is probably in frequencies, it's really hard to control and use them in manually calculating the frequency response; I suggest using the function filter; Look at the updated answer^ – Vidak Jun 3 '13 at 16:04

Vidak gives an excellent answer to your question #1, so I will just deal with #2. No, a high pass filter is not simply 1 - LPF. The reason why is that there is magnitude response and there is phase response. Although "1" and "LPF" do have the same magnitude response at the pass-band frequencies (1), they do not necessarily have the same phase response, so they could constructively or destructively interfere with each other.

The way to create a high-pass filter is to create a low-pass filter and then modulate it to the Nyquist frequency by multiplying it with a sinusoid whose frequency is the Nyquist frequency.

You can create a band-pass filter by the same method. Create a low-pass filter and then modulate it to the band-pass center frequency that you want.

• I don't understand the way to create the HPF. I have the LPF and then what to multiply with? I'm using the same M and Fs, if that what you meant. – jico Jun 2 '13 at 9:12

You already got two good answers, but I'd like to clarify the point about obtaining the impulse response of a highpass filter given the response of a lowpass filter. In your case (if you take the response suggested in Vidak's answer) the lowpass impulse response is symmetric w.r.t. its maximum value. This implies that the corresponding phase response is linear:

$$H(e^{j\theta})=A(\theta)e^{-j\theta (N-1)/2}\tag{1}$$

where $N$ is the (odd) filter length (i.e. the number of taps). In (1) it is assumed that the filter is causal and that the first non-zero filter tap is at $n=0$, i.e. the center of symmetry is at $n=(N-1)/2$ for odd $N$.

In this case you can indeed obtain the response of a highpass filter by a simple subtraction if you take the linear phase response into account:

$$h_{HP}(n)=\delta\left[ n-(N-1)/2\right ]-h_{LP}(n)$$

If you shift the lowpass filter impulse response to the left such that the center of symmetry is at the point $n=0$ (and the filter taps are non-zero between $n=-(N-1)/2$ and $n=(N-1)/2$) you obtain your original formula:

$$h_{HP}(n)=1-h_{LP}(n)$$

A formula of this type can be applied to all odd-length FIR filters with linear phase.

One of the simplest methods to build a low pass filter is using fir2 function in matlab. Here is the code which i use

fs=70MHz % Sampling freq = 70 MHz
fc=fs/(10);   % pass band corner frequency
fc=fs/(10);   % pass band corner frequency
fc1=fs/(8);   %stop band corner frequency
%change the scaling factor according to ur cutoff frequency
lpfLength=127; % Order/Number of Filter coefficients

% Low Pass Filterlpf=fir2(lpfLength,[0 fc fc1 fs/2]/(fs/2),[1 1 0 0]);

FFThn=fft(lpf,NFFT);    % FFT of filter h(n) i.e H(w). NFFT can be any value like 2048

magfhn=20*log10(abs(FFThn)).';    % Magnitude in dB
magfhn=circshift(magfhn,NFFT/2).';  %circshift is used if u want to plot it from -fs/2 to fs/2

figure;
plot((-NFFT/2:NFFT/2-1)*fs/NFFT,magfhn); %Freq. plot of LPF, If u want to plot from 0 - fs than use 0:NFFT
xlabel('Frequency (Hz)');
ylabel('Suppression (dB)');
title('Low Pass Prototype Filter')
axis([-fs/2 fs/2 -100 5]);

• "without builtin functions". i got my answer earlier anyway. – jico Jun 25 '13 at 8:59
fs=70MHz % Sampling freq = 70 MHz
fc=fs/(10);   % pass band corner frequency
fc=fs/(10);   % pass band corner frequency
fc1=fs/(8);   %stop band corner frequency
%change the scaling factor according to ur cutoff frequency
lpfLength=127; % Order/Number of Filter coefficients

% Low Pass Filterlpf=fir2(lpfLength,[0 fc fc1 fs/2]/(fs/2),[1 1 0 0]);

FFThn=fft(lpf,NFFT);    % FFT of filter h(n) i.e H(w). NFFT can be any value like 2048

magfhn=20*log10(abs(FFThn)).';    % Magnitude in dB
magfhn=circshift(magfhn,NFFT/2).';  %circshift is used if u want to plot it from -fs/2 to fs/2

figure;
plot((-NFFT/2:NFFT/2-1)*fs/NFFT,magfhn); %Freq. plot of LPF, If u want to plot from 0 - fs than use 0:NFFT
xlabel('Frequency (Hz)');
ylabel('Suppression (dB)');
title('Low Pass Prototype Filter')
axis([-fs/2 fs/2 -100 5]);

• a textual explanation of your code and/or more directly addressing the original questions would be helpful – PAK-9 Nov 20 '13 at 19:16