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

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:  1819 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:  1800 mal
filter_analyse.m
 Beschreibung:
Wird zur graphischen Darstellung der Filtereigenschaften benötigt.

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

Download
 Dateiname:  window_sinc_filter.m
 Dateigröße:  7.09 KB
 Heruntergeladen:  1731 mal
Private Nachricht senden Benutzer-Profile anzeigen


ckunkel

Gast


Beiträge: ---
Anmeldedatum: ---
Wohnort: ---
Version: ---
     Beitrag Verfasst am: 07.11.2018, 09:47     Titel: Signallänge
  Antworten mit Zitat      
Hallo!

Ich habe mir ein Minimalbeispiel der Filterung mit dem window-sinc_filter zusammengebaut und dabei festgestellt, dass das gefilterte Signal mehr Samples enthält als das Ausgangssignal. Das liegt wahrscheinlich an der hohen Filterordnung des Bandpasses vermute ich.
Das Problem der ganzen Geschichte liegt nun darin, dass es sich um ein Zeitsignal handelt, welches gefiltert werden soll. Wie kann ich nun also die Werte wieder richtig "sortieren" sodass ich im gefilterten Signal die korrekten Zeitpunkte erhalte?

Gruß
Christian

Code:

load sensorData
fs = Fs;

order = 500;
fc1 = 700;
fc2 = 800;
filter_type = 'bandpass';
window_type = '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
gefiltertes_Signal = FFT_Faltung(s1,h);
% äquivalent zur Faltung im Zeitbereich
gefiltertes_Signal2 = conv(s1,h);
% Faltung im Frequenzbereich ist jedoch bei höherer Filterordnung schneller als conv()

figure
hold on
plot(s1)
plot(gefiltertes_Signal)
 
 
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.