Hi,

Does any one know how specify the gain of a filter during design? It seems most design method does not provide unity gain. In fact, the always attenuates the amplitude of the input signal

Bayo

"Bayo " <ayb503@bham.ac.uk> wrote in message <hmo1ct$rjg$1@fred.mathworks.com>...
> Hi,
> 
> Does any one know how specify the gain of a filter during design? It seems most design method does not provide unity gain. In fact, the always attenuates the amplitude of the input signal
> 
> Bayo

Hi Bayo, presumably you are talking about unity gain in the passband. Can you give the MATLAB code you are using to generate your filter to demonstrate what gain you are getting in the passband? 

Wayne

"Wayne King" <wmkingty@gmail.com> wrote in message <hmo5e9$d5c$1@fred.mathworks.com>...
> "Bayo " <ayb503@bham.ac.uk> wrote in message <hmo1ct$rjg$1@fred.mathworks.com>...
> > Hi,
> > 
> > Does any one know how specify the gain of a filter during design? It seems most design method does not provide unity gain. In fact, the always attenuates the amplitude of the input signal
> > 
> > Bayo
> 
> Hi Bayo, presumably you are talking about unity gain in the passband. Can you give the MATLAB code you are using to generate your filter to demonstrate what gain you are getting in the passband? 
> 
> Wayne

Hi,

Thank you for your reply. 
The code below designs an elliptic highpass and lowpass filter concatenated to form a bandpass. I read in one of the help file that this concatenation work better than actual bandpass filter designs.

In addition, the output of this filter tails off at the end. Could you please enlighten me on how to fix this issue??

=====================================================
%Function: Design elliptic high pass and low pass filter using
%ellipord to generate the lowest no of filter order

function [LPf_eegData, HPf_eegData, BPf_eegData]=BPellipord(eegData)

%% Design LowPass Filter (17Hz)
Fs = 256; Rp=0.1; Rs=60; Wp = 16.5/(Fs/2); Ws = 17/(Fs/2); 
[n,Wp] = ellipord(Wp,Ws,Rp,Rs); %Generate minimum filter order for specs
[z,p,k] = ellip(n,Rp,Rs,Wp);
[sos,g] = zp2sos(z,p,k);	 %Convert to SOS form
Hd = dfilt.df2tsos(sos,g); %Create a dfilt object

%% Design Highpass filter (6 Hz)

Wn = 5.5/(Fs/2); Wm = 6/(Fs/2); %rising edge from 5 to 6
[n,Wm] = ellipord(Wn,Wm,Rp,Rs); %Generate minimum filter order for specs
[z,p,k] = ellip(n,Rp,Rs,Wm,'high');
[sos,g] = zp2sos(z,p,k);	 %Convert to SOS form
Hp = dfilt.df2tsos(sos,g); %Create a dfilt 
                  
 fvtool(Hp,Hd); 

%% Apply filter
[b,a] = sos2tf(Hd.sosMatrix,Hd.Scalevalues); 
[d,c] = sos2tf(Hp.sosMatrix,Hp.Scalevalues);
HPf_eegData = filtfilt(d,c,eegData);
LPf_eegData = filtfilt(b,a,eegData);
BPf_eegData = filtfilt(b,a,HPf_eegData);

==================================================

Thanks,

Bayo

"Bayo " <ayb503@bham.ac.uk> wrote in message <hmoffe$2hb$1@fred.mathworks.com>...
> "Wayne King" <wmkingty@gmail.com> wrote in message <hmo5e9$d5c$1@fred.mathworks.com>...
> > "Bayo " <ayb503@bham.ac.uk> wrote in message <hmo1ct$rjg$1@fred.mathworks.com>...
> > > Hi,
> > > 
> > > Does any one know how specify the gain of a filter during design? It seems most design method does not provide unity gain. In fact, the always attenuates the amplitude of the input signal
> > > 
> > > Bayo
> > 
> > Hi Bayo, presumably you are talking about unity gain in the passband. Can you give the MATLAB code you are using to generate your filter to demonstrate what gain you are getting in the passband? 
> > 
> > Wayne
> 
> Hi,
> 
> Thank you for your reply. 
> The code below designs an elliptic highpass and lowpass filter concatenated to form a bandpass. I read in one of the help file that this concatenation work better than actual bandpass filter designs.
> 
> In addition, the output of this filter tails off at the end. Could you please enlighten me on how to fix this issue??
> 
> =====================================================
> %Function: Design elliptic high pass and low pass filter using
> %ellipord to generate the lowest no of filter order
> 
> function [LPf_eegData, HPf_eegData, BPf_eegData]=BPellipord(eegData)
> 
> %% Design LowPass Filter (17Hz)
> Fs = 256; Rp=0.1; Rs=60; Wp = 16.5/(Fs/2); Ws = 17/(Fs/2); 
> [n,Wp] = ellipord(Wp,Ws,Rp,Rs); %Generate minimum filter order for specs
> [z,p,k] = ellip(n,Rp,Rs,Wp);
> [sos,g] = zp2sos(z,p,k);	 %Convert to SOS form
> Hd = dfilt.df2tsos(sos,g); %Create a dfilt object
> 
> %% Design Highpass filter (6 Hz)
> 
> Wn = 5.5/(Fs/2); Wm = 6/(Fs/2); %rising edge from 5 to 6
> [n,Wm] = ellipord(Wn,Wm,Rp,Rs); %Generate minimum filter order for specs
> [z,p,k] = ellip(n,Rp,Rs,Wm,'high');
> [sos,g] = zp2sos(z,p,k);	 %Convert to SOS form
> Hp = dfilt.df2tsos(sos,g); %Create a dfilt 
>                   
>  fvtool(Hp,Hd); 
> 
> %% Apply filter
> [b,a] = sos2tf(Hd.sosMatrix,Hd.Scalevalues); 
> [d,c] = sos2tf(Hp.sosMatrix,Hp.Scalevalues);
> HPf_eegData = filtfilt(d,c,eegData);
> LPf_eegData = filtfilt(b,a,eegData);
> BPf_eegData = filtfilt(b,a,HPf_eegData);
> 
> ==================================================
> 
> Thanks,
> 
> Bayo

Hi Bayo, 
I don't see anything obvously wrong here. Perhaps it has to do with the data. EEG data is usually nonstationary so that may be affecting things.  If I test your bandpass filter on some simulated stationary data, it works fine:

% your filter looks fine
hcas = dfilt.cascade(Hd,Hp);
fvtool(hcas)

n=1:1024;
% I know this is nothing like EEG data
eegData = cos(.1*pi*n)+randn(size(n));
HPf_eegData = filtfilt(d,c,eegData);
BPf_eegData = filtfilt(b,a,HPf_eegData);


Just out of curiosity if you design an FIR filter using fdesign.bandpass that matches your filter specs and use that in filtfilt() , do you get the same problem?

If you want you can send me an example EEG data file (along with the sampling frequency and what frequency band you'd like to pass) and I can look at it quickly. Not sure if I can say anything more substantive if I had the actual data, but...

Wayne

"Wayne King" <wmkingty@gmail.com> wrote in message <hmon19$m88$1@fred.mathworks.com>...
> "Bayo " <ayb503@bham.ac.uk> wrote in message <hmoffe$2hb$1@fred.mathworks.com>...
> > "Wayne King" <wmkingty@gmail.com> wrote in message <hmo5e9$d5c$1@fred.mathworks.com>...
> > > "Bayo " <ayb503@bham.ac.uk> wrote in message <hmo1ct$rjg$1@fred.mathworks.com>...
> > > > Hi,
> > > > 
> > > > Does any one know how specify the gain of a filter during design? It seems most design method does not provide unity gain. In fact, the always attenuates the amplitude of the input signal
> > > > 
> > > > Bayo
> > > 
> > > Hi Bayo, presumably you are talking about unity gain in the passband. Can you give the MATLAB code you are using to generate your filter to demonstrate what gain you are getting in the passband? 
> > > 
> > > Wayne
> > 
> > Hi,
> > 
> > Thank you for your reply. 
> > The code below designs an elliptic highpass and lowpass filter concatenated to form a bandpass. I read in one of the help file that this concatenation work better than actual bandpass filter designs.
> > 
> > In addition, the output of this filter tails off at the end. Could you please enlighten me on how to fix this issue??
> > 
> > =====================================================
> > %Function: Design elliptic high pass and low pass filter using
> > %ellipord to generate the lowest no of filter order
> > 
> > function [LPf_eegData, HPf_eegData, BPf_eegData]=BPellipord(eegData)
> > 
> > %% Design LowPass Filter (17Hz)
> > Fs = 256; Rp=0.1; Rs=60; Wp = 16.5/(Fs/2); Ws = 17/(Fs/2); 
> > [n,Wp] = ellipord(Wp,Ws,Rp,Rs); %Generate minimum filter order for specs
> > [z,p,k] = ellip(n,Rp,Rs,Wp);
> > [sos,g] = zp2sos(z,p,k);	 %Convert to SOS form
> > Hd = dfilt.df2tsos(sos,g); %Create a dfilt object
> > 
> > %% Design Highpass filter (6 Hz)
> > 
> > Wn = 5.5/(Fs/2); Wm = 6/(Fs/2); %rising edge from 5 to 6
> > [n,Wm] = ellipord(Wn,Wm,Rp,Rs); %Generate minimum filter order for specs
> > [z,p,k] = ellip(n,Rp,Rs,Wm,'high');
> > [sos,g] = zp2sos(z,p,k);	 %Convert to SOS form
> > Hp = dfilt.df2tsos(sos,g); %Create a dfilt 
> >                   
> >  fvtool(Hp,Hd); 
> > 
> > %% Apply filter
> > [b,a] = sos2tf(Hd.sosMatrix,Hd.Scalevalues); 
> > [d,c] = sos2tf(Hp.sosMatrix,Hp.Scalevalues);
> > HPf_eegData = filtfilt(d,c,eegData);
> > LPf_eegData = filtfilt(b,a,eegData);
> > BPf_eegData = filtfilt(b,a,HPf_eegData);
> > 
> > ==================================================
> > 
> > Thanks,
> > 
> > Bayo
> 
> Hi Bayo, 
> I don't see anything obvously wrong here. Perhaps it has to do with the data. EEG data is usually nonstationary so that may be affecting things.  If I test your bandpass filter on some simulated stationary data, it works fine:
> 
> % your filter looks fine
> hcas = dfilt.cascade(Hd,Hp);
> fvtool(hcas)
> 
> n=1:1024;
> % I know this is nothing like EEG data
> eegData = cos(.1*pi*n)+randn(size(n));
> HPf_eegData = filtfilt(d,c,eegData);
> BPf_eegData = filtfilt(b,a,HPf_eegData);
> 
> 
> Just out of curiosity if you design an FIR filter using fdesign.bandpass that matches your filter specs and use that in filtfilt() , do you get the same problem?
> 
> If you want you can send me an example EEG data file (along with the sampling frequency and what frequency band you'd like to pass) and I can look at it quickly. Not sure if I can say anything more substantive if I had the actual data, but...
> 
> Wayne

Hi Wayne,

Thanks for your reply.

I did design a FIR filter with the same specification but its performance was even worse. It literally attenuated most of the signal. 

I do not know how to attach files on here. However, I just sent some eeg data to personal email (wmkingty@gmail.com). Hope you do not mind.

Thank you .

Bayo

On 4 Mar, 11:18, "Bayo " <ayb...@bham.ac.uk> wrote:
> Hi,
>
> Does any one know how specify the gain of a filter during design? It seems most design method does not provide unity gain. In fact, the always attenuates the amplitude of the input signal

Cant's see what the problem might be, as that's what filters
are supposed to do...

Rune

If you want to gain your filter, just multiply the numerator of the transfer functions by the desired gain value. For example, to gain your low-pass filter by 30 dB you would multiply the b coeff by 10^(30/20). Try plotting the freqz of this:

g = 10^(30/20);
freqz(g.*b,a,2^15,Fs);

If I am understanding correctly, you believe the entire bandwidth of the signal is being attenuated. Are you looking at the results in the time domain only? If you look at the results in time domain only, it might look attenuated. But, if you analyze the results of the filter operation in the frequency domain, you will see the filter is working according to specifications. Add the following to the end of your script and run it again:

% Plot PSD
N = 2^10;
Pxx_eegData = pwelch(eegData,hann(N),0,N,Fs);
[Pxx_BPf_eegData f] = pwelch(BPf_eegData,hann(N),0,N,Fs);

figure;
hold on;
plot(f,10*log10(Pxx_eegData));
plot(f,10*log10(Pxx_BPf_eegData),'r:');
title('Welch Power Spectral Density Estimate');
xlabel('Frequency (Hz)');
ylabel('Power/frequency (dB/Hz)');
box on;
grid on;
legend('Orig EEG Data','BP Filtered EEG Data');
hold off;