이 질문을 팔로우합니다.
- 팔로우하는 게시물 피드에서 업데이트를 확인할 수 있습니다.
- 정보 수신 기본 설정에 따라 이메일을 받을 수 있습니다.
Signal amplitude decreases after FIR filter(windowing method)
조회 수: 13 (최근 30일)
이전 댓글 표시
Amy Lg
2022년 3월 10일
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
답변 (1개)
Star Strider
2022년 3월 10일
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
Amy Lg
2022년 3월 10일
편집: Amy Lg
2022년 3월 10일
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
2022년 3월 10일
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
2022년 3월 10일
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
2022년 3월 10일
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
2022년 3월 10일
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
2022년 3월 11일
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
2022년 3월 11일
편집: Amy Lg
2022년 3월 11일
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
2022년 3월 11일
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
2022년 3월 11일
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
2022년 3월 12일
What do you want to do with it?
Determine that, then design the filter around that requirement.
Star Strider
2022년 3월 12일
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
2022년 3월 12일
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
2022년 3월 12일
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
2022년 3월 12일
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
2022년 3월 12일
‘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
2022년 3월 12일
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
2022년 3월 12일
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
2022년 3월 19일
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
2022년 3월 19일
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
2022년 3월 19일
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
2022년 3월 19일
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.
참고 항목
카테고리
Help Center 및 File Exchange에서 Digital Filter Analysis에 대해 자세히 알아보기
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!오류 발생
페이지가 변경되었기 때문에 동작을 완료할 수 없습니다. 업데이트된 상태를 보려면 페이지를 다시 불러오십시오.
웹사이트 선택
번역된 콘텐츠를 보고 지역별 이벤트와 혜택을 살펴보려면 웹사이트를 선택하십시오. 현재 계신 지역에 따라 다음 웹사이트를 권장합니다:
또한 다음 목록에서 웹사이트를 선택하실 수도 있습니다.
사이트 성능 최적화 방법
최고의 사이트 성능을 위해 중국 사이트(중국어 또는 영어)를 선택하십시오. 현재 계신 지역에서는 다른 국가의 MathWorks 사이트 방문이 최적화되지 않았습니다.
미주
- América Latina (Español)
- Canada (English)
- United States (English)
유럽
- 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)
아시아 태평양
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)