You are now following this question
- You will see updates in your followed content feed.
- You may receive emails, depending on your communication preferences.
Signal amplitude decreases after FIR filter(windowing method)
12 views (last 30 days)
Show older comments
I am applying a low pass filter with the windowing method to my signal in order to truncate it in the frequency domain. The bandwidth of the filter should be 3GHz. I guess that works fine. The filter, however, seems to attenuate the amplitude of the signal. What do I need to do to get the same amplitude as the original signal after filtering?
I am using filters for the first time. I do appreciate any help or explanation.
sig = MY SIGNAL;
fs = 4000; % sampling freq. (GHz)
M = 400001; % signal length
% Filter parameters:
L = M; % filter length
fcut = 1.5; % cutoff frequency (GHz)
% Design the filter using the window method:
hsupp = (-(L-1)/2:(L-1)/2);
hideal = (2*fcut/fs)*sinc(2*fcut*hsupp/fs);
h = hamming(L)' .* hideal; % h is our filter
SIG_out = fft(sig); % signal
H = fft(h); % filter
FILT_OUT = SIG_out .* H;
filt_out = ifft(FILT_OUT);
relrmserr = norm(imag(filt_out))/norm(filt_out) % check... should be zero
Answers (1)
Star Strider
on 10 Mar 2022
I have no idea what the sampling frequency is, so I assume 10 GHz here, with a cutoff of 1.5 GHz. When I simulate it with the freqz function, it shows essentially no attenuation in the passband (as it should not), a stopband frequency of 1.5 GHz (as designed) and an appropriate transition region. Some parameters are not given (specifically ‘M’ and ‘fs’) so the filter behaviour here may be different from what your analysis shows. However on the basis of the freqz results, the problems you are seeing with it may simply be that your analysis — and not the filter design — is incorrect. (A larger value for ‘M’ will give a sharper cutoff and shorter transition region.)
M = 64; % Create 'M'
fs = 1E+10; % Create 'fs'
L = M; % filter length
fcut = 1.5E+9; % cutoff frequency (GHz)
% Design the filter using the window method:
hsupp = (-(L-1)/2:(L-1)/2);
hideal = (2*fcut/fs)*sinc(2*fcut*hsupp/fs);
h = hamming(L)' .* hideal; % h is our filter
figure
freqz(h, 1, 2^16, fs)
.
25 Comments
Amy Lg
on 10 Mar 2022
Edited: Amy Lg
on 10 Mar 2022
Thank you very much for your reply. The fs in my analysis is 4000GHz and M= 400001. The 'sig' is the output of some analysis and modeling in my system, and at this point, I need to truncate it in frequency domain based on 3GHz bandwidth. I do appreciate it if you could help me to figure out my mistake.
Star Strider
on 10 Mar 2022
My pleasure!
Actually, as designed it has a 1.5 Ghz bandwidth, and due to the filter length, a very short transition region. (I did not change either of those in this analysis.)
Changing those parameters —
M = 400001; % Use Supplied 'M'
fs = 4E+12; % Use Supplied 'fs'
L = M; % filter length
fcut = 1.5E+9; % cutoff frequency (GHz)
% Design the filter using the window method:
hsupp = (-(L-1)/2:(L-1)/2);
hideal = (2*fcut/fs)*sinc(2*fcut*hsupp/fs);
h = hamming(L)' .* hideal; % h is our filter
figure
freqz(h, 1, 2^16, fs)
set(subplot(2,1,1), 'XLim',[0 1E+10]) % Zoom Frequency Axis To Show Detail
set(subplot(2,1,2), 'XLim',[0 1E+10]) % Zoom Frequency Axis To Show Detail
That is a very long filter! Still, it demonstrates that it fulfills the design requirements.
If you are using the Signal Processing Toolbox functions, use the filtfilt function to do the actual filtering. The signal will need to be very long, or the filter length will need to be reduced (possibly significantly so). However it works efficiently as a much shorter filter, so that should not be a problem. If it is, then consider a IIR filter instead, most likely an elliptic design for numerical efficiency.
.
Amy Lg
on 10 Mar 2022
Thank you so much for your explanation. I am new to this topic and I wanna mke sure I am in the right direction. The all thing I need to do is truncate the signal in frequency domain. To this end, as you suggest, can I use IIR filters?
Star Strider
on 10 Mar 2022
My pleasure!
The filter will limit the frequency content of the filtered signal, so it will truncate it in the frequency domain, regardless of what filter design is used. Both the FIR and IIR filters will do this. The problem with the FIR filter of this length is that the signal being filtered must be more than twice the filter length. This is not the case for an IIR filter of any design (although they also have signal length requirements), however elliptic filters are computationally the most efficient. See the ellip function for details. (I always get the ‘z,p,k’ output, and then implement it with the zp2sos function and use that in the filtfilt call.)
.
Amy Lg
on 10 Mar 2022
Thank you so much. I am searcing to see how I should write the code for IIR filter. In the meantime, one more question. The output filtered signal shows oscillations in the bottom. Do you know why this is happening?
Star Strider
on 11 Mar 2022
It is difficult for me to figure out what is causing that, since I have no idea what the original signal is (so I could experiment with it) or how it was filtered (for example, using filtfilt).
It almost looks as though the filter is taking the derivative of the signal, and that is not something I would expect from a properly-designed lowpass filter. The filter is clearly ‘ringing’ in response to the initial spikes, however. An elliptic IIR filter might have a better response if implemented using the procedure I previously described.
Amy Lg
on 11 Mar 2022
Edited: Amy Lg
on 11 Mar 2022
I do appreciate your help and explanation. Based on your sugestion, I tried to implement the IIR filter as below. I was wondering if it is correct. Also, I am not sure based on what criteria we should choose appropriate values for filter order, passband ripple, and stopband ripple.
Regarding the ripple in passband, is it correct if I look at the abs of the 'filteredSignal' rather than 'filteredSignal'?
I would be so thankful if you could guide me in this regard.
n = ?; % filter order
fn = fs/2; % nyquist frequency
Rp = ?; % Passband Ripple (Attenuation)
Rs = ?; % Stopband Ripple (Attenuation)
fcut = 1.5; % cutoff frequency(GHz)
[z,p,k] = ellip(n,Rp,Rs,fcut/fn,'low'); % fcut has to be normalised to (divided by) the Nyquist frequency
[sos,g] = zp2sos(z,p,k);
filteredSignal = filtfilt(sos,g,sig_out);
Star Strider
on 11 Mar 2022
My pleasure!
You are close the the correct code, however it needs a few changes and additions.
I use the ellipord function to calculate the filter order, and generally do something like this to design it:
Fs = 1E+10; % Sampling Frequency
Fn = Fs/2; % Nyquist Frequency
Wp = 1.5E+9/Fn; % Passband Frequency (normalised)
Ws = 1.01*Wp; % Stopband Frequency (Normalised)
Rp = 1; % Passband Ripple (Attenuation)
Rs = 50; % Stopband Ripple (Attenuation)
[n,Wp] = ellipord(Wp,Ws,Rp,Rs); % Calculate Filter Order
[z,p,k] = ellip(n,Rp,Rs,Wp,'low'); % Zero-Pole-Gain Representation
[sos,g] = zp2sos(z,p,k); % Second-Order-Section Implementation
figure
freqz(sos, 2^16, Fs) % Filter Bode Plot
% filteredSignal = filtfilt(sos,g,sig_out); % Filter Signal
I commented the filtfilt call here so it would not execute when I run the rest of the code without the ‘sig_out’ vector. The syntax is correct.
Taking the abs fo ‘filteredSignal’ would simply force the negative parts to be positive. This would not desirable unless it is complex, and I doubt it would be complex. Even if it were complex, it would be better to consider the real and imaginary parts separately, rather than losing that information in combining them with abs.
Also, since the sampling frequency and filter frequency are in the GHz range, they must be specified accurately in order for the filter parameter matrix to be scaled correctly. Using 1.5 when the desired value is 1.5E+9 is not appropriate and will not produce the desired filter.
.
Star Strider
on 11 Mar 2022
I had forgotten what the sampling frequency was, and use the first value that I saw from my earlier post.
Use this filter instead —
Fs = 4E+12; % Sampling Frequency
Fn = Fs/2; % Nyquist Frequency
Wp = 1.5E+9/Fn; % Passband Frequency (normalised)
Ws = 1.01*Wp; % Stopband Frequency (Normalised)
Rp = 1; % Passband Ripple (Attenuation)
Rs = 50; % Stopband Ripple (Attenuation)
[n,Wp] = ellipord(Wp,Ws,Rp,Rs); % Calculate Filter Order
[z,p,k] = ellip(n,Rp,Rs,Wp,'low'); % Zero-Pole-Gain Representation
[sos,g] = zp2sos(z,p,k); % Second-Order-Section Implementation
figure
freqz(sos, 2^16, Fs) % Filter Bode Plot
set(subplot(2,1,1), 'XLim',[0 5E+9])
set(subplot(2,1,2), 'XLim',[0 5E+9])
% filteredSignal = filtfilt(sos,g,sig_out); % Filter Signal
Did you do a fft on the signal to determine what the frequency components are? If not, that needs to be done in order to design the filter correctly. The filter may have too narrow a passband, or the wrong passband, to filter your signal. I do not have your signal, so I cannot determine this.
.
Star Strider
on 12 Mar 2022
What do you want to do with it?
Determine that, then design the filter around that requirement.
Star Strider
on 12 Mar 2022
O.K. With that information, try this filter —
Fs = 4E+12; % Sampling Frequency
Fn = Fs/2; % Nyquist Frequency
Wp = 3E+9/Fn; % Passband Frequency (normalised)
Ws = 1.05*Wp; % Stopband Frequency (Normalised)
Rp = 1; % Passband Ripple (Attenuation)
Rs = 50; % Stopband Ripple (Attenuation)
[n,Wp] = ellipord(Wp,Ws,Rp,Rs); % Calculate Filter Order
[z,p,k] = ellip(n,Rp,Rs,Wp,'low'); % Zero-Pole-Gain Representation
[sos,g] = zp2sos(z,p,k); % Second-Order-Section Implementation
figure
freqz(sos, 2^16, Fs) % Filter Bode Plot
set(subplot(2,1,1), 'XLim',[0 5E+9])
set(subplot(2,1,2), 'XLim',[0 5E+9])
% filteredSignal = filtfilt(sos,g,sig_out); % Filter Signal
.
Amy Lg
on 12 Mar 2022
So you mean the cut-off frequency should be equal to the bandwidth limitation I have?
I though filter is centered around zero and we visualize just positive frequency? I was wrong?
Star Strider
on 12 Mar 2022
The filter eliminates everything (in a two-sided Fourier transform of the signal) less than -3 Ghz and greater than +3 GHz. (In a one-sided Fourier transform, everything greater than 3 GHz.)
I am not certain what ‘bandwidth limitation’ refers to here.
If the intent is to design a bandpass filter with a 3 GHz bandwidth, it will be necessary to know what the centre frequency is. The filter design after that requires only slight editing of my existing filter code.
Amy Lg
on 12 Mar 2022
The last designed filter is for a two-sided Fourier transform of the signal or one-sided Fourier transform?
The fourier transform of my signal is centered at zero. So, it has both positive and negative frequencies and its spectrum is two-sided. I got confused the designed filter with 'Wp = 3E+9/Fn' can be applied to my signal or not. The filter bandwidth should be 3GHz. Which one I should consider for wp, 3E+9/Fn or 1.5E+9/Fn?
Thank you
Star Strider
on 12 Mar 2022
‘The last designed filter is for a two-sided Fourier transform of the signal or one-sided Fourier transform?’
Actually, both. Using it on the signal and then taking the two-sided Fourier transform will demonstrate that there is energy only between -3 GHz and +3 GHz, with very little outside that region (since it will have been attenuated by 50 dB.). On a one-sided Fourier transform, it will only show energy between 0 and 3 GHz. The resulting time-domain signal is the same after filtering, the only difference is how the Fourier transform of it is plotted.
Again, I am not certain what ‘bandwidth limitation’ refers to here.
Amy Lg
on 12 Mar 2022
So, I am not able to use the latest designed filter. I should use Wp = 1.5E+9/Fn rather than Wp = 3E+9/Fn. I mean "bandwidth limitation" as follows:
I have a device with a 3GHz bandwidth limitation(ex. a photodetector). This must be applied to my signal. It's like a square-law demodulation. I should apply a LPF with 3GHz passband bandwidth to the intensity of my signal.
I really appreciate your patience for answering my questions.
Star Strider
on 12 Mar 2022
My pleasure!
I interpret that as:
Wp = 3E*9/Fn;
simply because ‘bandwidth’ is usually defined with respect to a single-sided Fourier transform.
So the intent is to design a lowpass filter with a 3 Ghz cutoff frequency.
Amy Lg
on 19 Mar 2022
I appreciate your help so much.
I still have problems with the result. Increasing my signal length to 4450001 and fs to 10000GHz results in an output similar to what is shown below. I do not know why somewhere the amplitude of the filtered signal increase and somewhere decrease. I think it is wrong to have oscillation in areas where the original signal is zero, because the low pass filter is supposed to smooth the signal and not adding new oscillation. I do appreciate if you could help me to correct the filter for my signal.
sig = MY SIGNAL; % Intensity of output signal(square-law device)
M = 4450001; % signal length
fs = 10000e9; % sampling freq. (10000GHz)
fcut = 3e9;
fn = fs/2; % Nyquist Frequency(Hz)
Wp = fcut/fn; % Passband Frequency (normalised)-Hz
Ws = 1.05*Wp; % Stopband Frequency (Normalised)-Hz
Rp = 1; % Passband Ripple (Attenuation)-dB
Rs = 50; % Stopband Ripple (Attenuation)-dB
[n,Wp] = ellipord(Wp,Ws,Rp,Rs); % Calculate Filter Order
[z,p,k] = ellip(n,Rp,Rs,Wp,'low'); % Zero-Pole-Gain Representation
[sos,g] = zp2sos(z,p,k); % Second-Order-Section Implementation
filteredSignal = filtfilt(sos,g,sig);
figure
freqz(sos, 2^16, fs) % Filter Bode Plot
Star Strider
on 19 Mar 2022
I have no idea what the problem is, since I do not understand what you want the filter to do.
I am not able to help you with this, so I will delete my answer in a few hours.
Amy Lg
on 19 Mar 2022
Why deleting the answer? Since so far I've learned a lot of new things with your help. I think it may help others who are new to this field like me. So, please do not delet your answer.
I think what I see in a square law device is the envelope of the signal (smoothed signal). So, I guess the filtered signal should follow the envelope of the original signal.
Star Strider
on 19 Mar 2022
O.K., I will keep it posted, at least for a while longer.
I simply have no idea how to solve the problems you have.
See Also
Categories
Find more on Single-Rate Filters in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!An Error Occurred
Unable to complete the action because of changes made to the page. Reload the page to see its updated state.
Select a Web Site
Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .
You can also select a web site from the following list
How to Get Best Site Performance
Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.
Americas
- América Latina (Español)
- Canada (English)
- United States (English)
Europe
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)
Asia Pacific
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)