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
% Berechnung der Filterkoeffizienten for i = 0:M
if2 * 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;
% 1. Filterstufe: Berechnung der Filterkoeffizienten for i = 0:M
if2 * 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
if2 * 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;
% 1. Filterstufe: Berechnung der Filterkoeffizienten for i = 0:M
if2 * 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
if2 * 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;
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
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?
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()
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
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.
RSS