Mein MATLAB Forum - goMatlab.de

Mein MATLAB Forum

 
Gast > Registrieren       Autologin?   
Bücher:

Fachkräfte:
Softwareentwickler (m/w) für automatische Codegenerierung
Softwareentwicklung mit MATLAB/Simulink und dSPACE TargetLink im Bereich Fahrwerkregelsysteme
Elektronische Fahrwerksysteme GmbH - Ingolstadt

Entwicklungsingenieur (m/w) für modellbasierte Softwareentwicklung
Modellbasierte Softwareentwicklung mit MATLAB/Simulink und dSPACE TargetLink
Elektronische Fahrwerksysteme GmbH - Ingolstadt

Test-Ingenieur (m/w) für Resimulation
Implementierung, Integration und Anpassung von sogenannten ADTF-Filtern
Automotive Safety Technologies GmbH - Ingolstadt

Resident-Ingenieur (w/m) Hardware-in-the-Loop-Simulation
Inbetriebnahme und Software-Anpassungen der HIL-Systeme
dSPACE GmbH - Wolfsburg

Entwicklungsingenieur (m/w) für Fahrdynamikregelsysteme
Entwicklung der für die Fahrdynamik relevanten Funktionen
Elektronische Fahrwerksysteme GmbH - Ingolstadt

weitere Angebote

Partner:


Vermarktungspartner


Forum
      Option
[Erweitert]
  • Diese Seite per Mail weiterempfehlen
     


Gehe zu:  
Neues Thema eröffnen Neue Antwort erstellen

Window-Sinc Filter

 

DSP
Forum-Meister

Forum-Meister



Beiträge: 2.117
Anmeldedatum: 28.02.11
Wohnort: ---
Version: R2014b
     Beitrag Verfasst am: 24.07.2011, 19:33     Titel: Window-Sinc Filter
  Antworten mit Zitat      
Diese Funktion berechnet die Koeffizienten eines Window-Sinc Filters. Es können insgesamt 4 verschiedene Filtertypen, Hochpass, Tiefpass, Bandpass und Bandsperre erstellt werden. Zu dem kann zwischen zwei Fenstern, Hamming und Blackman gewählt werden.

Code:

% Funktionsaufruf in Command Window oder m-file
order = 100; % Filterordnung -> 101 Koeffizienten
fc1 = 10; % 1. Grenzfrequenz
fc2 = 100; % 2. Grenzfrequenz
fs = 1000 % Abtastfrequenz in Hz
filter_type = 'bandpass';
window_typ = 'hamming';
analyse_plot = 'y'; % Filterfrequenzantwort darstellen
% Berechnung der Impulsantwort = Filterkoeffizienten
% Bandpass im Bereich 10...100 Hz
h = window_sinc_filter(order, fc1, fc2, fs, filter_type, window_type, analyse_plot)

% Signal mit dem Bandpass filtern
gefilteres Signal = FFT_Faltung(Signal,h) ;
% äquivalent zur Faltung im Zeitbereich
gefilteres Signal = conv(Signal,h)
% Faltung im Frequenzbereich ist jedoch bei höherer Filterordnung schneller als conv()
 


Zur Wahl des Fenstertyps:
Das Hamming Fenster weist eine 20% höhere Steilheit (vgl. Darstellung lineare Frequenzantwort) nach der Grenzfrequenz im Gegensatz zum Blackman Fenster auf. Das Blackman hat dagegen eine höhere Dämpfung im Sperrbereich und einen geringeren Rippel im Durchlassbereich 0.02% (Hamming: 0.2%).

Die Funktion zur Berechung des Fensters enthält noch weitere Fenstertypen...siehe:
http://www.gomatlab.de/viewtopic,p,71727.html#71727
Allerdings bieten sie keinen Vorteil im Bereich Dämpfung oder Steilheit gegenüber dem Hamming oder Blackman Fenster.

Die Grenzfrequenz/-en:
Die angegeben Grenzfrequenzen fc1 und fc2 werden nicht exakt erreicht. Die Abweichung hängt sowohl von der Filterordnung als auch von der Abtastfrequenz fs ab. Die tatsächlichen Grenzfrequenzen des Filters werden bei der Darstellung der Filteranalyse (anaylse_plot = 'y') bestimmt und ausgegeben.

Code:

function [output] = window_sinc_filter(order, fc1, fc2, fs, filter_type, window_type, analyse_plot)
% Die Funktion berechnet die Filterkoeffizienten eines Window-Sinc Filters.
% Hierzu werden die sinc Funktion sin(x)/x und ein wählbares Fenster verwendet.
%
% Inputparameter:
% order        = Filterordnung ->  Anzahl Filterkoeff. = order + 1
%                Bedingung: order = gerade ganze Zahl
%
% fc1 [Hz]     = cutoff frequency 1 = Grenzfrequenz der ersten Filterstufe
%                -> für Hochpass, Tiefpass, Bandpass und Bandsperre
%                Bedingung: 0 < fc1 < fs/2
%
% fc2 [Hz]     = cutoff frequency 2 = Grenzfrequenz der zweiten Filterstufe
%                -> neben fc1 zusätzlich nur für Bandpass und Bandsperre
%                Bedingung: 0 < fc1 < fc2 < fs/2
%
% fs [Hz]      = Sampling-frequency = Abtastfrequenz
%
% filter_type  = 'low' = Tiefpass, 'high' = Hochpass,, 'bandpass' = Bandpass
%                und 'bandreject' = Bandsperre
%
% window_type  = Fenstertyp: 'Hamming' oder 'Blackman'
%
% analyse_plot = Bodediagramm, Frequenzantwort linear und Sprungfunktion(nur Tiefpass)
%                des Filters darstellen: 'y' oder 'n'
%
% Outputparameter:
% output     = Filterkoeffizienten = Impulsantwort des Filters
%
% Quelle des Filter-Algorithmus:
% The Scientist and Engineer's Guide to Digital Signal Processing
% By Steven W. Smith, Ph.D.
% Chapter 16: Window-Sinc Filter
% http://www.dspguide.com/ch16.htm
% -------------------------------------------------------------------------

% Überprüfung der Inputparameter:
% Abtastfrequenz überprüfen
if (fs < 1) || (~isreal(fs))
    error(['Wrong sampling frequency fs [Hz]: ', sprintf('%8.1f',fs),' ...must be a real positiv number']);
end  
% Grenzfrequenz 1 überprüfen
Fc1 = fc1 / fs ; % Normierung auf fn (Fc1 = 0...0.5 fs)
if (Fc1 <= 0) || (Fc1 >= 0.5)
    error(['Wrong cutoff frequency fc1 [Hz]: ', sprintf('%8.1f',fc1),' ...choose: 0 < fc1 < ',sprintf('%8.1f',fs/2),' = fs/2']);
end
% Grenzfrequenz 2 überprüfen
% wird nur für Bandpass- und Bandsperrfilter benötigt
if (strcmp(filter_type, 'bandpass')) || (strcmp(filter_type, 'bandreject'))
    Fc2 = fc2 / fs ; % Normierung auf fn (Fc2 = 0...0.5 fs)
    if ((Fc2 <= 0) || (Fc2 >= 0.5) || (Fc2 <= Fc1))
        error(['Wrong cutoff frequency fc2 [Hz]: ', sprintf('%8.1f',fc2),' ...choose: ',sprintf('%8.1f',fc1) ' < fc2 < ',sprintf('%8.1f',fs/2),' = fs/2']);
    end
end
% Filterordnung überprüfen
M  = order;
if (mod(M, 2) ~= 0) % Filterordnung gerade?
    error(['Wrong Filterorder: ', sprintf('%8.1f',order),' ...must be even']);
end

% Window Typ überprüfen
if (~strcmp(window_type, 'hamming')) && (~strcmp(window_type, 'blackman'))
    error(['Unknown window type: ', window_type,' ...choose: hamming or blackman']);
end

% Auswahl Filtertyp
switch lower(filter_type)
   
    case 'high' % Hochpass - Filter
               
        output = zeros(M+1, 1); % Init
        % Fensterfunktion erstellen
        window = Fenster(M+1, window_type);
        % Berechnung der Filterkoeffizienten
        for i = 0:M
            if 2 * i == M
                B(i+1) = 2*pi*Fc1;
            else
                B(i+1) = sin(2*pi*Fc1 * (i-(M/2))) / (i-(M/2));
            end
            B(i+1) = B(i+1) * window(i+1);
        end                

        % Verstärkungsfaktor des Filters auf 1 normieren
        B = B./sum(B);
        % Filterkoeef. übergeben
        % Tiefpass in Hochpass durch Inversion des Spektrums wandeln
        output      = - B;
        output((M/2)+1) = output((M/2)+1) + 1;
       
    case 'low' % Tiefpass - Filter

        output = zeros(M+1, 1); % Init
        % Fensterfunktion erstellen
        window = Fenster(M+1, window_type);

        % Berechnung der Filterkoeffizienten
        for i = 0:M
            if 2 * i == M   % Multiplikation ist schneller als Division
                B(i+1) = 2*pi*Fc1;
            else
                B(i+1) = sin(2*pi*Fc1 * (i-(M/2))) / (i-(M/2));
            end
            B(i+1) = B(i+1) * window(i+1);
        end

        % Verstärkungsfaktor des Filters auf 1 normieren
        B      = B./sum(B);
        % Filterkoeef. übergeben
        output = B;

    case 'bandpass' % Bandpass - Filter

        output = zeros(M+1, 1); % Init
        % Fensterfunktion erstellen
        window = Fenster(M+1, window_type);
       
        % 1. Filterstufe: Berechnung der Filterkoeffizienten
        for i = 0:M
            if 2 * i == M   % Multiplikation ist schneller als Division
                A(i+1) = 2*pi*Fc1;
            else
                A(i+1) = sin(2*pi*Fc1 * (i-(M/2))) / (i-(M/2));
            end
            A(i+1) = A(i+1) * window(i+1);
        end

        % 2. Filterstufe: Berechnung der Filterkoeffizienten
        for i = 0:M
            if 2 * i == M
                B(i+1) = 2*pi*Fc2;
            else
                B(i+1) = sin(2*pi*Fc2 * (i-(M/2))) / (i-(M/2));
            end
            B(i+1) = B(i+1) * window(i+1);
        end

        % Verstärkungsfaktor des Filters auf 1 normieren
        A = A./sum(A);
        B = B./sum(B);
        % Tiefpass in Hochpass durch Inversion des Spektrums wandeln
        B          = - B;
        B((M/2)+1) = B((M/2)+1) + 1;
        output     = A + B;
        % Filterkoeef. übergeben
        % Bandsperre in Bandpass durch Inversion des Spektrums wandeln
        output          = - output;
        output((M/2)+1) = output((M/2)+1) + 1;

    case 'bandreject' % Bandsperr - Filter

        output = zeros(M+1, 1); % Init
        % Fensterfunktion erstellen
        window = Fenster(M+1, window_type);

        % 1. Filterstufe: Berechnung der Filterkoeffizienten
        for i = 0:M
            if 2 * i == M   % Multiplikation ist schneller als Division
                A(i+1) = 2*pi*Fc1;
            else
                A(i+1) = sin(2*pi*Fc1 * (i-(M/2))) / (i-(M/2));
            end
            A(i+1) = A(i+1) * window(i+1);
        end

        % 2. Filterstufe: Berechnung der Filterkoeffizienten
        for i = 0:M
            if 2 * i == M  
                B(i+1) = 2*pi*Fc2;
            else
                B(i+1) = sin(2*pi*Fc2 * (i-(M/2))) / (i-(M/2));
            end
            B(i+1) = B(i+1) * window(i+1);
        end

        % Verstärkungsfaktor des Filters auf 1 normieren
        A = A./sum(A);
        B = B./sum(B);
        % Tiefpass in Hochpass durch Inversion des Spektrums wandeln
        B          = - B;
        B((M/2)+1) = B((M/2)+1) + 1;
        % Filterkoeef. übergeben
        output = A + B;
       
    otherwise % ungültiger Filtertyp
        error(['Unknown filter type: ', filter_type ' ...choose: high, low, bandpass or bandreject']);
end % switch

% Filter Frequenzantwort darstellen?
% Eingabewert überprüfen
if (~strncmpi(analyse_plot, 'y', 1)) && (~strncmpi(analyse_plot, 'n', 1))
    error(['Unknown analyse_plot command: ', analyse_plot ,' ...choose: y or n']);
elseif (strncmpi(analyse_plot, 'y', 1))
    % Plots erstellen
     filter_analyse(output, fs, filter_type);
end
 


Die Algorithmus zur Berechung der Filterkoeffizienten stammt aus der folgenden Quelle:
The Scientist and Engineer's Guide to Digital Signal Processing
By Steven W. Smith, Ph.D.
Chapter 16: Window-Sinc Filter
http://www.dspguide.com/ch16.htm

Vielen Dank noch an Jan S für die Optimierung des Programmcodes!

Gruß,

DSP

Fenster.m
 Beschreibung:
Wird zur Erstellung des Fensters (Hamming oder Blackman) benötigt

Download
 Dateiname:  Fenster.m
 Dateigröße:  1.02 KB
 Heruntergeladen:  785 mal
FFT_Faltung.m
 Beschreibung:
Wird zur Filterung eines Signals mit dem erstellten Window-Sinc Filter benötigt

Download
 Dateiname:  FFT_Faltung.m
 Dateigröße:  1.02 KB
 Heruntergeladen:  844 mal
filter_analyse.m
 Beschreibung:
Wird zur graphischen Darstellung der Filtereigenschaften benötigt.

Download
 Dateiname:  filter_analyse.m
 Dateigröße:  6.76 KB
 Heruntergeladen:  794 mal
window_sinc_filter.m
 Beschreibung:
Funktion zur Berechnung der Filterkoeffizienten

Download
 Dateiname:  window_sinc_filter.m
 Dateigröße:  7.09 KB
 Heruntergeladen:  762 mal
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
.


goMatlab ist ein Teil des goForen-Labels
goForen.de goMATLAB.de goLaTeX.de goPCB.de


 Impressum  | Nutzungsbedingungen  | Datenschutz  | Werbung/Mediadaten | Studentenversion | FAQ | goMatlab RSS Button RSS


Copyright © 2007 - 2017 goMatlab.de | Dies ist keine offizielle Website der Firma The Mathworks
Partner: LabVIEWforum.de

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.