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

Question

related to Filtering voice signal from randomly added frequencies

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

Answer 1

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

Answer 2

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.

Answer 3

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.

Answer 4

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]); 

Answer 5

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]);