0
$\begingroup$

In an attempt to gain a better understanding of DSP, I want to create a very simple (1D) lowpass, or I think more correctly "stop-band" or "band-reject" filter, to filter out a single frequency, e.g., f2 in my code example below.

The idea is to build in in the frequency domain, by putting all ones values in every component of the filter's magnitude part EXCEPT of course in the position corresponding to f2. What i don't understand, is what phases am I would be supposed to put in the corresponding LP phase vector? (I think this is related to my problems of understanding phases in general)

For and ideal lowpass filter (brick-wall) I know the corresponding impulse response would be a sinc in the time domain. But e.g. if I built the magnitude part of the brickwall "manually" ones around 0 frequency and 0s for higher freq than the cutoff freq, I would have no idea what to put for the phase ?

My code:

clc;clear all;

fs=1000;
t=0:(1/fs):1-(1/fs);
f1=5;
f2=27;
A1=13;
A2=3;
foffset=3;
posFreqPos=f2+1;
negFreqPos=fs-f2+1;

s=A1*cos(2*pi*f1*t)+A2*cos(2*pi*f2*t);

S=fft(s);
Smag=abs(S);
Sphase=angle(S);

LP_mag=ones(size(S));
LP_mag(posFreqPos)=0;%filter out high freq f2
LP_mag(negFreqPos)=0;

%%%%%%%%%%%%%%%%%%%
% LPphase= ???
%%%%%%%%%%%%%%%%%%%

%the filter built from mag and phase
LPf=LP_mag.*exp(-i.*LPphase);

figure;
plot(LPf)

figure;
plot(s);
figure;
plot(Smag);

I have modified the code here, the first version didn't make sense sorry. It's really a question about the phase.

The idea is to build a filter such as: enter image description here

Following the convolution theorem which states that a point-wise multiplication in the frequency domain corresponds to a convolution in the time domain. Of course one simple solution which "roughly" works is to just set to zero the frequencies I want to remove (following the same idea as shown on the image ) in the magnitude image of the FFT of my image, the rebuild the image with

myImage_filtered = real(ifft(myImageMag_filtered.*exp(-i.*myImagephase))); %real(.) to avoid rounding error causing imaginary part to be non-zero

(image source:http://www.robots.ox.ac.uk/~az/lectures/ia/lect2.pdf )

|improve this question
$\endgroup$
  • 3
    $\begingroup$ That's a bad idea and won't work. See dsp.stackexchange.com/questions/6220/… $\endgroup$ – Hilmar Jul 30 at 18:16
  • 3
    $\begingroup$ Possible duplicate of Why is it a bad idea to filter by zeroing out FFT bins? $\endgroup$ – Marcus Müller Jul 30 at 21:40
  • $\begingroup$ @Hilmar yes sorry, I modified the code, the previous version didn't make much sense. I do this for understanding though, not for a "real" application. I would like to know which phase to put in the " %%% box" $\endgroup$ – Machupicchu Jul 31 at 7:59
  • $\begingroup$ @MarcusMüller I don't think so, since I am mostly interested about the phase, as I have said. $\endgroup$ – Machupicchu Jul 31 at 8:05
  • $\begingroup$ Sorry, your code is not doing anything particularly useful so the choice of phase will make no real difference to the outcome. You are asking the wrong question. Consider asking "what's the best way to remove a single frequency from a signal" and you will get answers that look very different from what you have. Filtering in the frequency domain is rather complicated and you need a good understanding of the math behind it before you can write code that works $\endgroup$ – Hilmar Jul 31 at 9:33
0
$\begingroup$

You should use the ‘symmetric’ flag in the ifft call. This will give you a real result. Or, you can do it yourself by imposing conjugate symmetry before taking the ifft (you’ll need to add another stopband centered at FS- f2). However, it’s still a bad way to implement a filter, as mentioned in the comments.

|improve this answer
$\endgroup$
  • $\begingroup$ I understand, sorry, i forgot about the negative frequencies. Please see my new piece of code. It's more about getting an intuitive understanding than building a "good" filter. What I am really asking is which phases should I put in the filter and why. $\endgroup$ – Machupicchu Jul 31 at 8:07

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.