WICHTIG: Der Betrieb von goMatlab.de wird privat finanziert fortgesetzt. - Mehr Infos...

Mein MATLAB Forum - goMatlab.de

Mein MATLAB Forum

  
Gast > Registrieren       Autologin?   

Partner:




Forum
      Option
[Erweitert]
  • Diese Seite per Mail weiterempfehlen
     


Gehe zu:  
Neues Thema eröffnen Neue Antwort erstellen

filtfilt() - Numerische Stabilität IIR Filter (Butterworth)
 

benwin
Forum-Anfänger

Forum-Anfänger

Beiträge: 11
Anmeldedatum: 28.01.13
Wohnort: ---
Version: R2011b
     Beitrag Verfasst am: 04.02.2013, 14:05     Titel: filtfilt() - Numerische Stabilität IIR Filter (Butterworth)
  Antworten mit Zitat      
Bin gerade dabei, mich in Matlab einzufuchsen, und habe damit angefangen, die verschieden Implementationen von Butterworth-Filtern zu vergleichen. Also Transfer Function-, Z-P-G- und fdesign-Varianten.
Der Code führt zu zwei recht anschaulichen Bildern (einfach mal durchlaufen lassen):
Code:
% --- INVESTIGATE DIFFERENTS FILTER IMPLEMENTATIONS ----------------------
% --- HERE: BUTTERWORTH FILTERS ONLY -------------------------------------

% OBJECT ORIENTED (DSP System Toolbox)
% use actual frequencies
% fdesign: user specifies the desired characteristics of filter and
%          saves those specifications in a Filter Specification Object
% design: user implements filter as dfilt object, then used to filter data

% NOT OBJECT ORIENTED (Signal Processing Toolbox)
% use normalized frequencies
% butter(), cheby1/2(), ellip(), equiripple() ... (all IIR filters)

% Convention: unit frequency = Nyquist frequency = sample frequency / 2
%             1 = fn = fs / 2
% Example:  fs=1000Hz, fn=500Hz, fcut=29Hz, fcut_unit=29/500=0.058
% ------------------------------------------------------------------------
scrsz = get(0,'ScreenSize');

fs = 1000;          % sampling frequency
flop = 29;          % low pass frequency
fhip = 2.5;         % high pass frequency
fb = [fhip flop];   % band pass frequencies vector

fn = fs/2;          % Nyquist frequency
flopn = flop/fn;    % normalized low pass frequency
fhipn = fhip/fn;    % normalized high pass frequency
fbn = fb./fn;       % normalized band pass frequencies vector


% ------------------------------------------------------------------------
% FILTER DESIGN
n = 6;       % order of Butterworth filter
% butter():  returns an order 2*n digital bandpass filter if
%            fbn is a two-element vector
% fdesign(): n determines both the numerator and denominator orders
%            for IIR filters

% Bandpass Filter -------------------------------
fiType = 'bandpass';
% butter () - Transfer Function design - numerically unstable
[bandB,bandA] = butter(n,fbn,fiType);
fiBandTF = dfilt.df2(bandB,bandA);
% butter () - Zero-Pole-Gain design - numerically stable
[z,p,k] = butter(n,fbn,fiType);
[sos,g] = zp2sos(z,p,k);
fiBandZPG = dfilt.df2sos(sos,g);
% fdesign(), design() - normalized frequencies
fiSpec = fdesign.bandpass('N,F3dB1,F3dB2',2*n,fhipn,flopn);
fiBandDn = design(fiSpec,'butter');

% Lowpass Filter --------------------------------
fiType = 'low';
% butter () - Transfer Function design - numerically unstable
[lowB,lowA] = butter(n,flopn,fiType);
fiLowTF = dfilt.df2(lowB,lowA);
% butter () - Zero-Pole-Gain design - numerically stable
[z,p,k] = butter(n,flopn,fiType);
[sos,g] = zp2sos(z,p,k);
fiLowZPG = dfilt.df2sos(sos,g);
% fdesign(), design() - normalized frequencies
fiSpec = fdesign.lowpass('N,F3dB',n,flopn);
fiLowDn = design(fiSpec,'butter');

% Highpass Filter --------------------------------
fiType = 'high';
% butter () - Transfer Function design - numerically unstable
[highB,highA] = butter(n,fhipn,fiType);
fiHighTF = dfilt.df2(highB,highA);
% butter () - Zero-Pole-Gain design - numerically stable
[z,p,k] = butter(n,fhipn,fiType);
[sos,g] = zp2sos(z,p,k);
fiHighZPG = dfilt.df2sos(sos,g);
% fdesign(), design() - normalized frequencies
fiSpec = fdesign.highpass('N,F3dB',n,fhipn);
fiHighDn = design(fiSpec,'butter');




% ------------------------------------------------------------------------
% TIME SIGNAL FILTERING
rng default;
duration = 0.3;
ns = duration*fs;               % number of samples
t = linspace(0,duration,ns);    % time vector
y = ecg(ns)'+0.25*randn(ns,1);  % electro-cardiogram + noise

% Bandpass Filter -------------------------------
% phase distorted
yBandTF = filter(fiBandTF,y);
yBandTF1D = filter(bandB,bandA,y);
yBandZPG = filter(fiBandZPG,y);
yBandDn = filter(fiBandDn,y);
% zero-phase
yBandTF0p = filtfilt(fiBandTF.Numerator,1,y);
yBandTF1D0p = filtfilt(bandB,bandA,y);
yBandZPG0p = filtfilt(fiBandZPG.sosMatrix,fiBandZPG.ScaleValues,y);
yBandDn0p = filtfilt(fiBandDn.sosMatrix,fiBandDn.ScaleValues,y);

% Lowpass Filter --------------------------------
% phase distorted
yLowTF = filter(fiLowTF,y);
yLowTF1D = filter(lowB,lowA,y);
yLowZPG = filter(fiLowZPG,y);
yLowDn = filter(fiLowDn,y);
% zero-phase
yLowTF0p = filtfilt(fiLowTF.Numerator,1,y);
yLowTF1D0p = filtfilt(lowB,lowA,y);
yLowZPG0p = filtfilt(fiLowZPG.sosMatrix,fiLowZPG.ScaleValues,y);
yLowDn0p = filtfilt(fiLowDn.sosMatrix,fiLowDn.ScaleValues,y);

% Highpass Filter -------------------------------
% phase distorted
yHighTF = filter(fiHighTF,y);
yHighTF1D = filter(highB,highA,y);
yHighZPG = filter(fiHighZPG,y);
yHighDn = filter(fiHighDn,y);
% zero-phase
yHighTF0p = filtfilt(fiHighTF.Numerator,1,y);
yHighTF1D0p = filtfilt(highB,highA,y);
yHighZPG0p = filtfilt(fiHighZPG.sosMatrix,fiHighZPG.ScaleValues,y);
yHighDn0p = filtfilt(fiHighDn.sosMatrix,fiHighDn.ScaleValues,y);




% ------------------------------------------------------------------------
% TIME SIGNAL PLOTS
minT = 0; maxT = 0.3; % time
minY = -4; maxY = 4; % values
titPX = minT; titPY = 4.02; % title plot
titRX = -0.05; titRY = 0; % title row

setBorderX = @() xlim([minT,maxT]);
setBorderY = @() ylim([minY,maxY]);
titProps = {'HorizontalAlignment','left','VerticalAlignment','bottom','Fontsize',8,'Position',[titPX titPY]};
rowTxt = @(s) text(titRX,titRY,s,'Rotation',90,'HorizontalAlignment','center');

figSig = figure('Name','Filtered Time Signals', ...
    'Position',[scrsz(3)*0.05 scrsz(4)*0.1 scrsz(3)*0.8 scrsz(4)*0.8]);
% Band-Passed Signals ---------------------------
subplot(4,3,1); plot(t,y,'k'); hold on;
plot(t,yBandTF,t,yBandTF0p,'linewidth',1.5);
setBorderX(); setBorderY();
title({'filter(fiBandTF,x)','filtfilt(fiBandTF.Numerator,1,x)'},titProps{:});
rowTxt('TF butter()');

subplot(4,3,4); plot(t,y,'k'); hold on;
plot(t,yBandTF1D,t,yBandTF1D0p,'linewidth',1.5);
setBorderX(); setBorderY();
title({'filter(lowB,lowA,x)','filtfilt(bandB,bandA,x)'},titProps{:});
rowTxt('TF butter()');

subplot(4,3,7);plot(t,y,'k');hold on;
plot(t,yBandZPG,t,yBandZPG0p,'linewidth',1.5);
setBorderX(); setBorderY();
title({'filter(fiBandZPG,x)','filtfilt(fiBandZPG.sosMatrix,fiBandZPG.ScaleValues,x)'},titProps{:});
rowTxt('ZPG butter()');

subplot(4,3,10);plot(t,y,'k');hold on;
plot(t,yBandDn,t,yBandDn0p,'linewidth',1.5);
setBorderX(); setBorderY();
title({'filter(fiBandDn,x)','filtfilt(fiBandDn.sosMatrix,fiBandDn.ScaleValues,x)'},titProps{:});
rowTxt('fdesign()');

% Low-Passed Signals ----------------------------
subplot(4,3,2);plot(t,y,'k');hold on;
plot(t,yLowTF,t,yLowTF0p,'linewidth',1.5);
title({'filter(fiLowTF,x)','filtfilt(fiLowTF.Numerator,1,x)'},titProps{:});
setBorderX(); setBorderY();

subplot(4,3,5);plot(t,y,'k');hold on;
plot(t,yLowTF1D,t,yLowTF1D0p,'linewidth',1.5);
title({'filter(lowB,lowA,x)','filtfilt(lowB,lowA,x)'},titProps{:});
setBorderX(); setBorderY();

subplot(4,3,8);plot(t,y,'k');hold on;
plot(t,yLowZPG,t,yLowZPG0p,'linewidth',1.5);
title({'filter(fiLowZPG,x)','filtfilt(fiLowZPG.sosMatrix,fiLowZPG.ScaleValues,x)'},titProps{:});
setBorderX(); setBorderY();

subplot(4,3,11);plot(t,y,'k');hold on;
plot(t,yLowDn,t,yLowDn0p,'linewidth',1.5);
title({'filter(fiLowDn,x)','filtfilt(fiLowDn.sosMatrix,fiLowDn.ScaleValues,x)'},titProps{:});
setBorderX(); setBorderY();

% High-Passed Signals ----------------------------
subplot(4,3,3);plot(t,y,'k');hold on;
plot(t,yHighTF,t,yHighTF0p,'linewidth',1.5);
title({'filter(fiHighTF,x)','filtfilt(fiHighTF.Numerator,1,x)'},titProps{:});
setBorderX(); setBorderY();

subplot(4,3,6);plot(t,y,'k');hold on;
plot(t,yHighTF1D,t,yHighTF1D0p,'linewidth',1.5);
title({'filter(highB,highA,x)','filtfilt(highB,highA,x)'},titProps{:});
setBorderX(); setBorderY();

subplot(4,3,9);plot(t,y,'k');hold on;
plot(t,yHighZPG,t,yHighZPG0p,'linewidth',1.5);
title({'filter(fiHighZPG,x)','filtfilt(fiHighZPG.sosMatrix,fiHighZPG.ScaleValues,x)'},titProps{:});
setBorderX(); setBorderY();

subplot(4,3,12);plot(t,y,'k');hold on;
plot(t,yHighDn,t,yHighDn0p,'linewidth',1.5);
title({'filter(fiHighDn,x)','filtfilt(fiHighDn.sosMatrix,fiHighDn.ScaleValues,x)'},titProps{:});
setBorderX(); setBorderY();




% ------------------------------------------------------------------------
% FILTER PLOTS
minR = 0; maxR = 2; % response
minF = 1e-3; maxF = 1; % frequency
titPX = minF; titPY = 2.01; % title plot
titRX = 2e-4; titRY = 1; % title row
stabX = 1.5*minF; stabYr1 = 0.9*maxR; stabYr2 = 0.78*maxR;

setBorderX = @() xlim([minF,maxF]);
setBorderY = @() ylim([minR,maxR]);
titProps = {'HorizontalAlignment','left','VerticalAlignment','bottom','Fontsize',8,'Position',[titPX titPY]};
rowTxt = @(s) text(titRX,titRY,s,'Rotation',90,'HorizontalAlignment','center');
stabTxt = @(f,posY) text(stabX,posY,['isstable(' inputname(1) '): ' int2str(isstable(f))],'HorizontalAlignment','left','Fontsize',8);
stabTxt2 = @(fB,fA,posY) text(stabX,posY,['isstable(' inputname(1) ',' inputname(2) '): ' ...
    int2str(isstable(dfilt.df2(fB,fA)))],'HorizontalAlignment','left','Fontsize',8);

figFil = figure('Name','Filters', ...
    'Position',[scrsz(3)*0.05 scrsz(4)*0.1 scrsz(3)*0.8 scrsz(4)*0.8]);
% Bandpass Filter -------------------------------
[magBandTF,fBandTF]=freqz(fiBandTF,5000);
subplot(4,3,1);semilogx(fBandTF,abs(magBandTF),'b');hold on;
title({'freqz(fiBandTF)'},titProps{:});
setBorderX(); setBorderY();
stabTxt(fiBandTF,stabYr1);
rowTxt('TF butter()');

[magBandTF1D,fBandTF1D]=freqz(bandB,bandA,5000);
subplot(4,3,4);semilogx(fBandTF1D,abs(magBandTF1D),'b');hold on;
title({'freqz(bandB,bandA)'},titProps{:});
setBorderX(); setBorderY();
stabTxt2(bandB,bandA,stabYr1);
rowTxt('TF butter()');

[magBandZPG,fBandZPG]=freqz(fiBandZPG,5000); % regular design
[bandZPG_B,bandZPG_A]= sos2tf(fiBandZPG.sosMatrix,fiBandZPG.ScaleValues);
[magBandZPG_TF,fBandZPG_TF]= freqz(bandZPG_B,bandZPG_A,5000); % transformed desing for filtfilt()
subplot(4,3,7);semilogx(fBandZPG,abs(magBandZPG),fBandZPG_TF,abs(magBandZPG_TF));hold on;
title({'freqz(fiBandZPG)','freqz(bandZPG\_B,bandZPG\_A)'},titProps{:});
setBorderX(); setBorderY();
stabTxt(fiBandZPG,stabYr1);
stabTxt2(bandZPG_B,bandZPG_A,stabYr2);
rowTxt('ZPG butter()');

[magBandDn,fBandDn]=freqz(fiBandDn,5000); % regular design
[bandDn_B,bandDn_A]= sos2tf(fiBandDn.sosMatrix,fiBandDn.ScaleValues);
[magBandDn_TF,fBandDn_TF]= freqz(bandDn_B,bandDn_A,5000); % transformed desing for filtfilt()
subplot(4,3,10);semilogx(fBandDn,abs(magBandDn),fBandDn_TF,abs(magBandDn_TF));hold on;
title({'freqz(fiBandDn)','freqz(bandDn\_B,bandDn\_A)'},titProps{:});
setBorderX(); setBorderY();
stabTxt(fiBandDn,stabYr1);
stabTxt2(bandDn_B,bandDn_A,stabYr2);
rowTxt('fdesign()');

% Lowpass Filter --------------------------------
[magLowTF,fLowTF]=freqz(fiLowTF,5000);
subplot(4,3,2);semilogx(fLowTF,abs(magLowTF),'b');hold on;
title({'freqz(fiLowTF)'},titProps{:});
setBorderX(); setBorderY();
stabTxt(fiLowTF,stabYr1);

[magLowTF1D,fLowTF1D]=freqz(lowB,lowA,5000);
subplot(4,3,5);semilogx(fLowTF1D,abs(magLowTF1D),'b');hold on;
title({'freqz(lowB,lowA)'},titProps{:});
setBorderX(); setBorderY();
stabTxt2(lowB,lowA,stabYr1);

[magLowZPG,fLowZPG]=freqz(fiLowZPG,5000); % regular design
[lowZPG_B,lowZPG_A]= sos2tf(fiLowZPG.sosMatrix,fiLowZPG.ScaleValues);
[magLowZPG_TF,fLowZPG_TF]= freqz(lowZPG_B,lowZPG_A,5000); % transformed desing for filtfilt()
subplot(4,3,8);semilogx(fLowZPG,abs(magLowZPG),fLowZPG_TF,abs(magLowZPG_TF));hold on;
title({'freqz(fiLowZPG)','freqz(lowZPG\_B,lowZPG\_A)'},titProps{:});
setBorderX(); setBorderY();
stabTxt(fiLowZPG,stabYr1);
stabTxt2(lowZPG_B,lowZPG_A,stabYr2);

[magLowDn,fLowDn]=freqz(fiLowDn,5000); % regular design
[lowDn_B,lowDn_A]= sos2tf(fiLowDn.sosMatrix,fiLowDn.ScaleValues);
[magLowDn_TF,fLowDn_TF]= freqz(lowDn_B,lowDn_A,5000); % transformed desing for filtfilt()
subplot(4,3,11);semilogx(fLowDn,abs(magLowDn),fLowDn_TF,abs(magLowDn_TF));hold on;
title({'freqz(fiLowDn)','freqz(lowDn\_B,lowDn\_A)'},titProps{:});
setBorderX(); setBorderY();
stabTxt(fiLowDn,stabYr1);
stabTxt2(lowDn_B,lowDn_A,stabYr2);

% Highpass Filter -------------------------------
[magHighTF,fHighTF]=freqz(fiHighTF,5000);
subplot(4,3,3);semilogx(fHighTF,abs(magHighTF),'b');hold on;
title({'freqz(fiHighTF)'},titProps{:});
setBorderX(); setBorderY();
stabTxt(fiHighTF,stabYr1);

[magHighTF1D,fHighTF1D]=freqz(highB,highA,5000);
subplot(4,3,6);semilogx(fHighTF1D,abs(magHighTF1D),'b');hold on;
title({'freqz(highB,highA)'},titProps{:});
setBorderX(); setBorderY();
stabTxt2(highB,highA,stabYr1);

[magHighZPG,fHighZPG]=freqz(fiHighZPG,5000); % regular design
[highZPG_B,highZPG_A]= sos2tf(fiHighZPG.sosMatrix,fiHighZPG.ScaleValues);
[magHighZPG_TF,fHighZPG_TF]= freqz(highZPG_B,highZPG_A,5000); % transformed desing for filtfilt()
subplot(4,3,9);semilogx(fHighZPG,abs(magHighZPG),fHighZPG_TF,abs(magHighZPG_TF));hold on;
title({'freqz(fiHighZPG)','freqz(highZPG\_B,highZPG\_A)'},titProps{:});
setBorderX(); setBorderY();
stabTxt(fiHighZPG,stabYr1);
stabTxt2(highZPG_B,highZPG_A,stabYr2);

[magHighDn,fHighDn]=freqz(fiHighDn,5000); % regular design
[highDn_B,highDn_A]= sos2tf(fiHighZPG.sosMatrix,fiHighZPG.ScaleValues);
[magHighDn_TF,fHighDn_TF]= freqz(highDn_B,highDn_A,5000); % transformed desing for filtfilt()
subplot(4,3,12);semilogx(fHighDn,abs(magHighDn),fHighDn_TF,abs(magHighDn_TF));hold on;
title({'freqz(fiHighDn)','freqz(highZPG\_B,highZPG\_A)'},titProps{:});
setBorderX(); setBorderY();
stabTxt(fiHighDn,stabYr1);
stabTxt2(highDn_B,highDn_A,stabYr2);

Jetzt zur Frage:
Dass Highpass-Filter mit niedriger Cutoff-Frequenz numerisch instabil sein können ist klar. Bandpass-Filter entsprechend. Die Transfer-Function-Form ist dabei natürlich noch instabiler als die Alternativen - was blöd ist, wenn man das Signal ohne Phasenversatz filtern möchte und auf filtfilt() angewiesen ist. filtfilt() akzeptiert ja leider keine andere Form.

Wie kann man erfahren, ob ein Filter stabil unter filtfilt() ist?
Wenn ich beispielsweise umwandle:
Code:
[highZPG_B,highZPG_A]= sos2tf(fiHighZPG.sosMatrix,fiHighZPG.ScaleValues);
ergeben
Code:
isstable(fiHighZPG)
isstable(highZPG_B,highZPG_A)
beide 1, obwohl
Code:
y=filtfilt(highZPG_B,highZPG_A,x)
klar instabil ist (siehe rechte Spalte im Zeitsignalplot).

Woher kommt diese Instabilität wenn nicht aus den Polen?

Grüße,
Benjamin
Private Nachricht senden Benutzer-Profile anzeigen


Jan S
Moderator

Moderator

Beiträge: 11.057
Anmeldedatum: 08.07.10
Wohnort: Heidelberg
Version: 2009a, 2016b
     Beitrag Verfasst am: 04.02.2013, 20:44     Titel: Re: filtfilt() - Numerische Stabilität IIR Filter (Butterwo
  Antworten mit Zitat      
Hallo Bejamin,

Beim Versuch Deinen Code laufen zu lassen, scheitert mein Matlab am Fehlern der Funktion RNG. Zur eigentlichen Frage fehlt mir also leider das Anschauungsmaterial.

Der Unterschied von FILTFILT zum doppelt angewendeten FILTER ist lediglich die besondere Methode, die versucht an den Rändern das Einschwingverhalten zu reduzieren. Es bleibt aber natürlich sowieso fraglich, ob automatisch erzeugte zusätzliche Signal-Stücke eine realistische Verlängerung des echten Signals darstellen, oder ob dies Artefakte sind. Es gibt also durchaus auch Fälle, in denen FILTFILT wissenschaftlich sinnvoll hierdurch ersetzt werden kann:
Code:
filtered = fliplr(filter(fliplr(filter(signal))));

Gruß, Jan
Private Nachricht senden Benutzer-Profile anzeigen
 
Harald
Forum-Meister

Forum-Meister

Beiträge: 24.492
Anmeldedatum: 26.03.09
Wohnort: Nähe München
Version: ab 2017b
     Beitrag Verfasst am: 04.02.2013, 20:55     Titel:
  Antworten mit Zitat      
Hallo Jan,

rng wurde in 2011a (glaube ich) eingeführt, um das Kontrollieren des Zufallszahlengenerators zu erleichtern.
http://www.mathworks.com/help/relea.....om-number-generation.html

Zum Testen sollte man den rng-Befehl hier auch weglassen können.

Grüße,
Harald
Private Nachricht senden Benutzer-Profile anzeigen
 
benwin
Themenstarter

 Forum-Anfänger

Forum-Anfänger

Beiträge: 11
Anmeldedatum: 28.01.13
Wohnort: ---
Version: R2011b
     Beitrag Verfasst am: 06.02.2013, 09:19     Titel:
  Antworten mit Zitat      
Der Hinweis auf die Besonderheit von filtfilt() war sehr hilfreich. Ich hab jetzt zusätzlich noch die flipud()-Variante mit in mein Beispiel aufgenommen (falls es jemanden interessiert ... und ja, die rng-Zeile kann auskommentiert werden.)
Anscheinend ist es wirklich so, dass das Verhalten, dass ich als numerische Instabiliät im Ergebnis von filtfilt() interpretiert habe, tatsächlich nur ein Artefakt und mit den Ergebnisse von isstable() konsistent ist.

Danke nochmal!
Code:
% --- INVESTIGATE DIFFERENTS FILTER IMPLEMENTATIONS ----------------------
% --- HERE: BUTTERWORTH FILTERS ONLY -------------------------------------

% OBJECT ORIENTED (DSP System Toolbox)
% use actual frequencies
% fdesign: user specifies the desired characteristics of filter and
%          saves those specifications in a Filter Specification Object
% design: user implements filter as dfilt object, then used to filter data

% NOT OBJECT ORIENTED (Signal Processing Toolbox)
% use normalized frequencies
% butter(), cheby1/2(), ellip(), equiripple() ... (all IIR filters)

% Convention: unit frequency = Nyquist frequency = sample frequency / 2
%             1 = fn = fs / 2
% Example:  fs=1000Hz, fn=500Hz, fcut=29Hz, fcut_unit=29/500=0.058
% ------------------------------------------------------------------------
scrsz = get(0,'ScreenSize');

fs = 1000;          % sampling frequency
flop = 29;          % low pass frequency
fhip = 2.5;         % high pass frequency
fb = [fhip flop];   % band pass frequencies vector

fn = fs/2;          % Nyquist frequency
flopn = flop/fn;    % normalized low pass frequency
fhipn = fhip/fn;    % normalized high pass frequency
fbn = fb./fn;       % normalized band pass frequencies vector


% ------------------------------------------------------------------------
% FILTER DESIGN
n = 6;       % order of Butterworth filter
% butter():  returns an order 2*n digital bandpass filter if
%            fbn is a two-element vector
% fdesign(): n determines both the numerator and denominator orders
%            for IIR filters

% Bandpass Filter -------------------------------
fiType = 'bandpass';
% butter () - Transfer Function design
[bandB,bandA] = butter(n,fbn,fiType);
fiBandTF = dfilt.df2(bandB,bandA);
% butter () - Zero-Pole-Gain design
[z,p,k] = butter(n,fbn,fiType);
[sos,g] = zp2sos(z,p,k);
fiBandZPG = dfilt.df2sos(sos,g);
% fdesign(), design() - here: normalized frequencies
fiSpec = fdesign.bandpass('N,F3dB1,F3dB2',2*n,fhipn,flopn);
fiBandDn = design(fiSpec,'butter');

% Lowpass Filter --------------------------------
fiType = 'low';
% butter () - Transfer Function design
[lowB,lowA] = butter(n,flopn,fiType);
fiLowTF = dfilt.df2(lowB,lowA);
% butter () - Zero-Pole-Gain design
[z,p,k] = butter(n,flopn,fiType);
[sos,g] = zp2sos(z,p,k);
fiLowZPG = dfilt.df2sos(sos,g);
% fdesign(), design() - here: normalized frequencies
fiSpec = fdesign.lowpass('N,F3dB',n,flopn);
fiLowDn = design(fiSpec,'butter');

% Highpass Filter --------------------------------
fiType = 'high';
% butter () - Transfer Function design
[highB,highA] = butter(n,fhipn,fiType);
fiHighTF = dfilt.df2(highB,highA);
% butter () - Zero-Pole-Gain design
[z,p,k] = butter(n,fhipn,fiType);
[sos,g] = zp2sos(z,p,k);
fiHighZPG = dfilt.df2sos(sos,g);
% fdesign(), design() - here: normalized frequencies
fiSpec = fdesign.highpass('N,F3dB',n,fhipn);
fiHighDn = design(fiSpec,'butter');




% ------------------------------------------------------------------------
% TIME SIGNAL FILTERING
rng default;
duration = 0.3;
ns = duration*fs;               % number of samples
t = linspace(0,duration,ns);    % time vector
y = ecg(ns)'+0.25*randn(ns,1);  % electro-cardiogram + noise

% Bandpass Filter -------------------------------
% phase distorted
yBandTF = filter(fiBandTF,y);
yBandTF1D = filter(bandB,bandA,y);
yBandZPG = filter(fiBandZPG,y);
yBandDn = filter(fiBandDn,y);
% zero-phase, filtfilt()
yBandTF0p = filtfilt(fiBandTF.Numerator,1,y);
yBandTF1D0p = filtfilt(bandB,bandA,y);
yBandZPG0p = filtfilt(fiBandZPG.sosMatrix,fiBandZPG.ScaleValues,y);
yBandDn0p = filtfilt(fiBandDn.sosMatrix,fiBandDn.ScaleValues,y);
% zero-phase, flipud(filter(flipud(filter())))
yBandTF0p2 = flipud(filter(fiBandTF,flipud(yBandTF)));
yBandTF1D0p2 = flipud(filter(bandB,bandA,flipud(yBandTF1D)));
yBandZPG0p2 = flipud(filter(fiBandZPG,flipud(yBandZPG)));
yBandDn0p2 = flipud(filter(fiBandDn,flipud(yBandDn)));

% Lowpass Filter --------------------------------
% phase distorted
yLowTF = filter(fiLowTF,y);
yLowTF1D = filter(lowB,lowA,y);
yLowZPG = filter(fiLowZPG,y);
yLowDn = filter(fiLowDn,y);
% zero-phase, filtfilt()
yLowTF0p = filtfilt(fiLowTF.Numerator,1,y);
yLowTF1D0p = filtfilt(lowB,lowA,y);
yLowZPG0p = filtfilt(fiLowZPG.sosMatrix,fiLowZPG.ScaleValues,y);
yLowDn0p = filtfilt(fiLowDn.sosMatrix,fiLowDn.ScaleValues,y);
% zero-phase, flipud(filter(flipud(filter())))
yLowTF0p2 = flipud(filter(fiLowTF,flipud(yLowTF)));
yLowTF1D0p2 = flipud(filter(lowB,lowA,flipud(yLowTF1D)));
yLowZPG0p2 = flipud(filter(fiLowZPG,flipud(yLowZPG)));
yLowDn0p2 = flipud(filter(fiLowDn,flipud(yLowDn)));

% Highpass Filter -------------------------------
% phase distorted
yHighTF = filter(fiHighTF,y);
yHighTF1D = filter(highB,highA,y);
yHighZPG = filter(fiHighZPG,y);
yHighDn = filter(fiHighDn,y);
% zero-phase, filtfilt()
yHighTF0p = filtfilt(fiHighTF.Numerator,1,y);
yHighTF1D0p = filtfilt(highB,highA,y);
yHighZPG0p = filtfilt(fiHighZPG.sosMatrix,fiHighZPG.ScaleValues,y);
yHighDn0p = filtfilt(fiHighDn.sosMatrix,fiHighDn.ScaleValues,y);
% zero-phase, flipud(filter(flipud(filter())))
yHighTF0p2 = flipud(filter(fiHighTF,flipud(yHighTF)));
yHighTF1D0p2 = flipud(filter(highB,highA,flipud(yHighTF1D)));
yHighZPG0p2 = flipud(filter(fiHighZPG,flipud(yHighZPG)));
yHighDn0p2 = flipud(filter(fiHighDn,flipud(yHighDn)));




% ------------------------------------------------------------------------
% TIME SIGNAL PLOTS
minT = 0; maxT = 0.3; % time
minY = -4; maxY = 5; % values
titPX = 0.005; titPY = maxY-0.01*(maxY-minY); % title plot
titRX = -0.05; titRY = 0; % title row

ploBorderX = @() xlim([minT,maxT]);
ploBorderY = @() ylim([minY,maxY]);
ploFontSize = @() set(gca,'FontSize',8);
ploSize = @() set(gca,'Position',[1 1 1.1 1.1].*get(gca,'Position'));
ploFormat = @() cellfun(@feval,{ploBorderX ploBorderY ploFontSize ploSize});
titProps = {'HorizontalAlignment','left','VerticalAlignment','top','Position',[titPX titPY]};
rowTxt = @(s) text(titRX,titRY,s,'Rotation',90,'HorizontalAlignment','center');
legDefault = @() legend('#1','#2','#3');
legFormat = @() set(legend(gca),'Position',[0 0.93 0.1 0.07],'Box','off','Color','none');
legPlot = @() cellfun(@feval,{legDefault legFormat});

figSig = figure('Name','Filtered Time Signals', ...
    'Position',[scrsz(3)*0.05 scrsz(4)*0.1 scrsz(3)*0.8 scrsz(4)*0.8]);
% Band-Passed Signals ---------------------------
subplot(4,3,1);
plot(t,yBandTF,t,yBandTF0p,t,yBandTF0p2,'linewidth',1.5); hold on;
plot(t,y,'k');
ploFormat();
title({'#1: filter(fiBandTF,x)','#2: filtfilt(fiBandTF.Numerator,1,x)','#3: flipud(filter(fiBandTF,flipud(#1)))'},titProps{:});
rowTxt('TF butter()');
legPlot();

subplot(4,3,4);
plot(t,yBandTF1D,t,yBandTF1D0p,t,yBandTF1D0p2,'linewidth',1.5); hold on;
plot(t,y,'k');
ploFormat();
title({'#1: filter(lowB,lowA,x)','#2: filtfilt(bandB,bandA,x)','#3: flipud(filter(lowB,lowA,flipud(#1)))'},titProps{:});
rowTxt('TF butter()');

subplot(4,3,7);
plot(t,yBandZPG,t,yBandZPG0p,t,yBandZPG0p2,'linewidth',1.5);hold on;
plot(t,y,'k');
ploFormat();
title({'#1: filter(fiBandZPG,x)','#2: filtfilt(fiBandZPG.sosMatrix,fiBandZPG.ScaleValues,x)','#3: flipud(filter(fiBandZPG,flipud(#1)))'},titProps{:});
rowTxt('ZPG butter()');

subplot(4,3,10);
plot(t,yBandDn,t,yBandDn0p,t,yBandDn0p2,'linewidth',1.5);hold on;
plot(t,y,'k');
ploFormat();
title({'#1: filter(fiBandDn,x)','#2: filtfilt(fiBandDn.sosMatrix,fiBandDn.ScaleValues,x)','#3: flipud(filter(fiBandDn,flipud(#1)))'},titProps{:});
rowTxt('fdesign()');

% Low-Passed Signals ----------------------------
subplot(4,3,2);
plot(t,yLowTF,t,yLowTF0p,t,yLowTF0p2,'linewidth',1.5);hold on;
plot(t,y,'k');
ploFormat();
title({'#1: filter(fiLowTF,x)','#2: filtfilt(fiLowTF.Numerator,1,x)','#3: flipud(filter(fiLowTF,flipud(#1)))'},titProps{:});

subplot(4,3,5);
plot(t,yLowTF1D,t,yLowTF1D0p,t,yLowTF1D0p2,'linewidth',1.5);hold on;
plot(t,y,'k');
ploFormat();
title({'#1: filter(lowB,lowA,x)','#2: filtfilt(lowB,lowA,x)','#3: flipud(filter(lowB,lowA,flipud(#1)))'},titProps{:});

subplot(4,3,8);
plot(t,yLowZPG,t,yLowZPG0p,t,yLowZPG0p2,'linewidth',1.5);hold on;
plot(t,y,'k');
ploFormat();
title({'#1: filter(fiLowZPG,x)','#2: filtfilt(fiLowZPG.sosMatrix,fiLowZPG.ScaleValues,x)','#3: flipud(filter(fiLowZPG,flipud(#1)))'},titProps{:});

subplot(4,3,11);
plot(t,yLowDn,t,yLowDn0p,t,yLowDn0p2,'linewidth',1.5);hold on;
plot(t,y,'k');
ploFormat();
title({'#1: filter(fiLowDn,x)','#2: filtfilt(fiLowDn.sosMatrix,fiLowDn.ScaleValues,x)','#3: flipud(filter(fiLowDn,flipud(#1)))'},titProps{:});

% High-Passed Signals ----------------------------
subplot(4,3,3);
plot(t,yHighTF,t,yHighTF0p,t,yHighTF0p2,'linewidth',1.5);hold on;
plot(t,y,'k');
ploFormat();
title({'#1: filter(fiHighTF,x)','#2: filtfilt(fiHighTF.Numerator,1,x)','#3: flipud(filter(fiHighTF,flipud(#1)))'},titProps{:});

subplot(4,3,6);
plot(t,yHighTF1D,t,yHighTF1D0p,t,yHighTF1D0p2,'linewidth',1.5);hold on;
plot(t,y,'k');
ploFormat();
title({'#1: filter(highB,highA,x)','#2: filtfilt(highB,highA,x)','#3: flipud(filter(highB,highA,flipud(#1)))'},titProps{:});

subplot(4,3,9);
plot(t,yHighZPG,t,yHighZPG0p,t,yHighZPG0p2,'linewidth',1.5);hold on;
plot(t,y,'k');
ploFormat();
title({'#1: filter(fiHighZPG,x)','#2: filtfilt(fiHighZPG.sosMatrix,fiHighZPG.ScaleValues,x)','#3: flipud(filter(fiHighZPG,flipud(#1)))'},titProps{:});

subplot(4,3,12);plot(t,y,'k');hold on;
plot(t,yHighDn,t,yHighDn0p,t,yHighDn0p2,'linewidth',1.5);hold on;
plot(t,y,'k');
ploFormat();
title({'#1: filter(fiHighDn,x)','#2: filtfilt(fiHighDn.sosMatrix,fiHighDn.ScaleValues,x)','#3: flipud(filter(fiHighDn,flipud(#1)))'},titProps{:});




% ------------------------------------------------------------------------
% FILTER PLOTS
minR = 0; maxR = 2; % response
minF = 1e-3; maxF = 1; % frequency
titPX = 1.2*minF; titPY = maxR-0.01*(maxR-minR); % title plot
titRX = 2e-4; titRY = 1; % title row
stabX = 1.2*minF; stabYr1 = 0.5*maxR; stabYr2 = 0.38*maxR; stabY = 0.01*maxR;

ploBorderX = @() xlim([minF,maxF]);
ploBorderY = @() ylim([minR,maxR]);
ploFontSize = @() set(gca,'FontSize',8);
ploSize = @() set(gca,'Position',[1 1 1.1 1.1].*get(gca,'Position'));
ploFormat = @() cellfun(@feval,{ploBorderX ploBorderY ploFontSize ploSize});
titProps = {'HorizontalAlignment','left','VerticalAlignment','top','Position',[titPX titPY]};
rowTxt = @(s) text(titRX,titRY,s,'Rotation',90,'HorizontalAlignment','center');
legDefault = @() legend('#1','#2');
legFormat = @() set(legend(gca),'Position',[0 0.95 0.1 0.05],'Box','off','Color','none');
legPlot = @() cellfun(@feval,{legDefault legFormat});
stabiTxt = @(i) text(stabX,stabY,['isstable(#1): ' int2str(i)],'HorizontalAlignment','left','VerticalAlignment','bottom','Fontsize',8);
stabiTxt2 = @(i1,i2) text(stabX,stabY,{['isstable(#1): ' int2str(i1)],['isstable(#2): ' int2str(i2)]},'HorizontalAlignment','left','VerticalAlignment','bottom','Fontsize',8);

figFil = figure('Name','Filters', ...
    'Position',[scrsz(3)*0.05 scrsz(4)*0.1 scrsz(3)*0.8 scrsz(4)*0.8]);
% Bandpass Filter -------------------------------
[magBandTF,fBandTF]=freqz(fiBandTF,5000);
subplot(4,3,1);
semilogx(fBandTF,abs(magBandTF));hold on;
ploFormat();
title({'#1: freqz(fiBandTF)'},titProps{:});
stabiTxt(isstable(fiBandTF));
rowTxt('TF butter()');

[magBandTF1D,fBandTF1D]=freqz(bandB,bandA,5000);
subplot(4,3,4);
semilogx(fBandTF1D,abs(magBandTF1D));hold on;
ploFormat();
title({'#1: freqz(bandB,bandA)'},titProps{:});
stabiTxt(isstable(dfilt.df2(bandB,bandA)));
rowTxt('TF butter()');

[magBandZPG,fBandZPG]=freqz(fiBandZPG,5000); % regular design
[bandZPG_B,bandZPG_A]= sos2tf(fiBandZPG.sosMatrix,fiBandZPG.ScaleValues);
[magBandZPG_TF,fBandZPG_TF]= freqz(bandZPG_B,bandZPG_A,5000); % transformed design for filtfilt()
subplot(4,3,7);
semilogx(fBandZPG,abs(magBandZPG),fBandZPG_TF,abs(magBandZPG_TF));hold on;
ploFormat();
title({'#1: freqz(fiBandZPG)','#2: freqz(bandZPG\_B,bandZPG\_A)'},titProps{:});
stabiTxt2(isstable(fiBandZPG),isstable(dfilt.df2(bandZPG_B,bandZPG_A)));
rowTxt('ZPG butter()');
legPlot();

[magBandDn,fBandDn]=freqz(fiBandDn,5000); % regular design
[bandDn_B,bandDn_A]= sos2tf(fiBandDn.sosMatrix,fiBandDn.ScaleValues);
[magBandDn_TF,fBandDn_TF]= freqz(bandDn_B,bandDn_A,5000); % transformed design for filtfilt()
subplot(4,3,10);
semilogx(fBandDn,abs(magBandDn),fBandDn_TF,abs(magBandDn_TF));hold on;
ploFormat();
title({'#1: freqz(fiBandDn)','#2: freqz(bandDn\_B,bandDn\_A)'},titProps{:});
stabiTxt2(isstable(fiBandDn),isstable(dfilt.df2(bandDn_B,bandDn_A)));
rowTxt('fdesign()');

% Lowpass Filter --------------------------------
[magLowTF,fLowTF]=freqz(fiLowTF,5000);
subplot(4,3,2);
semilogx(fLowTF,abs(magLowTF),'b');hold on;
ploFormat();
title({'#1: freqz(fiLowTF)'},titProps{:});
stabiTxt(isstable(fiLowTF));

[magLowTF1D,fLowTF1D]=freqz(lowB,lowA,5000);
subplot(4,3,5);
semilogx(fLowTF1D,abs(magLowTF1D),'b');hold on;
ploFormat();
title({'#1: freqz(lowB,lowA)'},titProps{:});
stabiTxt(isstable(dfilt.df2(lowB,lowA)));

[magLowZPG,fLowZPG]=freqz(fiLowZPG,5000); % regular design
[lowZPG_B,lowZPG_A]= sos2tf(fiLowZPG.sosMatrix,fiLowZPG.ScaleValues);
[magLowZPG_TF,fLowZPG_TF]= freqz(lowZPG_B,lowZPG_A,5000); % transformed design for filtfilt()
subplot(4,3,8);
semilogx(fLowZPG,abs(magLowZPG),fLowZPG_TF,abs(magLowZPG_TF));hold on;
ploFormat();
title({'#1: freqz(fiLowZPG)','#2: freqz(lowZPG\_B,lowZPG\_A)'},titProps{:});
stabiTxt2(isstable(fiLowZPG),isstable(dfilt.df2(lowZPG_B,lowZPG_A)));

[magLowDn,fLowDn]=freqz(fiLowDn,5000); % regular design
[lowDn_B,lowDn_A]= sos2tf(fiLowDn.sosMatrix,fiLowDn.ScaleValues);
[magLowDn_TF,fLowDn_TF]= freqz(lowDn_B,lowDn_A,5000); % transformed design for filtfilt()
subplot(4,3,11);
semilogx(fLowDn,abs(magLowDn),fLowDn_TF,abs(magLowDn_TF));hold on;
ploFormat();
title({'#1: freqz(fiLowDn)','#2: freqz(lowDn\_B,lowDn\_A)'},titProps{:});
stabiTxt2(isstable(fiLowDn),isstable(dfilt.df2(lowDn_B,lowDn_A)));

% Highpass Filter -------------------------------
[magHighTF,fHighTF]=freqz(fiHighTF,5000);
subplot(4,3,3);
semilogx(fHighTF,abs(magHighTF),'b');hold on;
ploFormat();
title({'#1: freqz(fiHighTF)'},titProps{:});
stabiTxt(isstable(fiHighTF));

[magHighTF1D,fHighTF1D]=freqz(highB,highA,5000);
subplot(4,3,6);
semilogx(fHighTF1D,abs(magHighTF1D));hold on;
ploFormat();
title({'#1: freqz(highB,highA)'},titProps{:});
stabiTxt(isstable(dfilt.df2(highB,highA)));

[magHighZPG,fHighZPG]=freqz(fiHighZPG,5000); % regular design
[highZPG_B,highZPG_A]= sos2tf(fiHighZPG.sosMatrix,fiHighZPG.ScaleValues);
[magHighZPG_TF,fHighZPG_TF]= freqz(highZPG_B,highZPG_A,5000); % transformed design for filtfilt()
subplot(4,3,9);
semilogx(fHighZPG,abs(magHighZPG),fHighZPG_TF,abs(magHighZPG_TF));hold on;
ploFormat();
title({'#1: freqz(fiHighZPG)','#2: freqz(highZPG\_B,highZPG\_A)'},titProps{:});
stabiTxt2(isstable(fiHighZPG),isstable(dfilt.df2(highZPG_B,highZPG_A)));

[magHighDn,fHighDn]=freqz(fiHighDn,5000); % regular design
[highDn_B,highDn_A]= sos2tf(fiHighZPG.sosMatrix,fiHighZPG.ScaleValues);
[magHighDn_TF,fHighDn_TF]= freqz(highDn_B,highDn_A,5000); % transformed design for filtfilt()
subplot(4,3,12);
semilogx(fHighDn,abs(magHighDn),fHighDn_TF,abs(magHighDn_TF));hold on;
ploFormat();
title({'#1: freqz(fiHighDn)','#2: freqz(highZPG\_B,highZPG\_A)'},titProps{:});
stabiTxt2(isstable(fiHighDn),isstable(dfilt.df2(highDn_B,highDn_A)));
Private Nachricht senden Benutzer-Profile anzeigen
 
Neues Thema eröffnen Neue Antwort erstellen



Einstellungen und Berechtigungen
Beiträge der letzten Zeit anzeigen:

 Du kannst Beiträge in dieses Forum schreiben.
Du kannst auf Beiträge in diesem Forum antworten.
Du kannst deine Beiträge in diesem Forum nicht bearbeiten.
Du kannst deine Beiträge in diesem Forum nicht löschen.
Du kannst an Umfragen in diesem Forum nicht mitmachen.
Du kannst Dateien in diesem Forum posten
Du kannst Dateien in diesem Forum herunterladen
.





 Impressum  | Nutzungsbedingungen  | Datenschutz | FAQ | goMatlab RSS Button RSS

Hosted by:


Copyright © 2007 - 2024 goMatlab.de | Dies ist keine offizielle Website der Firma The Mathworks

MATLAB, Simulink, Stateflow, Handle Graphics, Real-Time Workshop, SimBiology, SimHydraulics, SimEvents, and xPC TargetBox are registered trademarks and The MathWorks, the L-shaped membrane logo, and Embedded MATLAB are trademarks of The MathWorks, Inc.