Φίλτα σε Matlab

[gI@]

Γεια σας και καλή χρονιά, θα μπορούσε κάποιον να μου πει πως μπορώ να κάνω ένα αναλογικό Bandstop Butterworth φίλτρο στο matlab. Έχω γράψει κάποιο κώδικα αλλά ουδέν αποτέλεσμα.

Το σήμα μου έχει τις εξής συχνότητες (108Ηz 144z 180hz) και θέλω να αφαιρέσω την μεσαία, πιθανών κάνω λάθος στις συχνότητες Wp kai Ws...

Όποια βοήθεια θα ήταν καλή,
Σας ευχαριστώ

[osek]

Αν δοκιμάσεις να γράψεις στο Command Window, help buttord, θα δεις τα ορίσματά του. Το Wp και το Ws, είναι κυκλικές συχνότητες και όχι συχνότητες... Αν θες, δείξε τον κώδικά σου.

[gI@]

Ναι το ξέρω, το έχω κάνει αυτό
ο κώδικα μου είναι ο εξης:
Rp = 1;
Rs = 60;
Wp = 2*pi*[120,170];
Ws = 2*pi*[140,150];


[N,Wn]=buttord(Wp,Ws, Rp,Rs,'s')
Fc=Wn/(2*pi)
[num,den]=butter(N,Wn,'stop','s');

Δεν ξέρω εάν έχω επιλέξει σωστά όρια, το πρόβλημα εμφανίζετε αφού πάω να το εφαρμόσω στο σήμα με την "lsim"

[osek]

Θα πρότεινα αρχικά το φίλτρο να γίνει κάπως έτσι:

Rp = 1;
Rs = 60;
Wp = 2*pi*[140, 150];
Ws = 2*pi*[120, 170];
[N,Wn] = buttord(Wp,Ws, Rp,Rs,'s');
[num,den] = butter(N,Wn,'stop','s');
[H,W] = freqs(num, den);
plot(W/(2*pi), abs(H));

Και μετά να μεταφέρεις το σήμα σου στο πεδίο τον συχνοτήτων με φουριέ.

[gI@]

osek σε ευχαριστώ για το ενδιαφέρον σου αλλά ούτε αυτό δούλεψε, καταρχήν στο
[H,W] = freqs(num, den);plot(W/(2*pi), abs(H));η απόκριση του φίλτρου ξεκινάει από το 1 αντί για 0 και μου δίνει θετικά πράγμα που ενισχύει όλες τις συχνότητες στην έξοδοΣου παραθέτω όλο των κώδικα για να δεις τυχών λάθη μου :fs=400
Rp = 1;
Rs = 60;
Wp = 2*pi*[140, 150];
Ws = 2*pi*[120, 170];
[N,Wn] = buttord(Wp,Ws, Rp,Rs,'s');
[num,den] = butter(N,Wn,'stop','s');
[H,W] = freqs(num, den);
plot(W/(2*pi), abs(H));
sys=tf(num,den);
z=lsim(sys,x,t);
plot(t,z);
n=4096;
F=fft(z,n);
F0 = fftshift(F);
pedio = (-n/2:n/2-1)*(fs/n);
E0 = F0.*conj(F0)/n;
plot(pedio,E0)

[osek]

καταρχήν στο [H,W] = freqs(num, den);plot(W/(2*pi), abs(H));η απόκριση του φίλτρου ξεκινάει από το 1 αντί για 0 και μου δίνει θετικά πράγμα που ενισχύει όλες τις συχνότητες στην έξοδο

Φίλε Δημήτρη,
κάπου μπερδεφτήκαμε... η απόκριση ενός bandstop φίλτρου είναι κάπως έτσι,
27573
ξεκινάει από μέγιστο, περνάει στο ελάχιστο (περιοχή αποκοπής) και ξαναγυρνάει στο μέγιστο (περιοχή διάβασης). Αυτό που θέλεις κάπως έτσι δεν πρέπει να είναι;

Ο παρακάτω κώδικας βγάζει το εξής διάγραμμα.

Rp = 1;
Rs = 60;
Wp = 2*pi*[140, 150];
Ws = 2*pi*[120, 170];
[N,Wn] = buttord(Wp,Ws, Rp,Rs,'s');
[num,den] = butter(N,Wn,'stop','s');
[H,W] = freqs(num, den);
plot(W/(2*pi), abs(H));

27572

[gI@]

Ναι το ξέρω αυτό, sorry εγώ έκανα λάθος για το από πού θα αρχίζει.
Και εγώ αυτό βγάζω, αλλά όταν εφαρμόζετε στο σήμα μου δεν κάνει καμία αποκοπή στην μεσαία συχνότητα.

[gI@]

και επίσης θέλω να σε ρωτήσω για το κενό που εμφανίζετε όταν κάνεις ζουμ στο μηδέν. Γιατί συμβαίνει αυτό;

[osek]

και επίσης θέλω να σε ρωτήσω για το κενό που εμφανίζετε όταν κάνεις ζουμ στο μηδέν. Γιατί συμβαίνει αυτό;
Τώρα που το είπες το παρατήρησα και εγώ. Αλλά δεν ξέρω γιατί συμβαίνει αυτό.


Δες τον παρακάτω κώδικα με τα αντίστοιχα σχόλια. Ελπίζω να βοηθήσει.


clc

Rp = 1;
Rs = 60;
Wp = 2*pi*[140, 150];
Ws = 2*pi*[120, 170];
[N,Wn] = buttord(Wp,Ws, Rp,Rs,'s');
[num,den] = butter(N,Wn,'stop','s');
[H,W] = freqs(num, den);
plot(W/(2*pi), abs(H));
%Apokrisi filtrou


fs = 10*280; %fs=10*fmax
n=0:1/fs:1; %vima deigmatolipsias
x=sin(2*pi*44*n) + sin(2*pi*145*n) + sin(2*pi*280*n); %sima
D = 512; %o FFT tha ginei me 512 deigmata
X = abs(fft(x,D)); %vriskoume ton FFT tou x
f = ((0:D-1)/D*fs); %kanonikopoioume os pros ti suxnotita
figure;plot(f,X);
%fasmatiki apokrisi simatos


sys=tf(num,den);
z=lsim(sys,x,n); %pername to sima apo to filtro
D = 512; %o FFT tha ginei me 512 deigmata
Z = abs(fft(z,D)); %vriskoume ton FFT tis eksodou tou filtrou
f = ((0:D-1)/D*fs); %kanonikopoioume os pros ti suxnotita
figure;plot(f,Z);
%fasmatiki apokrisi simatos, meta to filtro

[gI@]

Σε ευχαριστώ πολύ τώρα δούλεψε μια χαρά και βρήκα ποιο είναι το πρόβλημα που είχα,

Η fs θέλει να είναι αρκετά ποιο μεγάλη, από την Nyquist συχνότητα. Την έβαλα 2000 και δούλεψε μια χαρά.

Και πάλι ευχαριστώ, δεν το πιστεύω ότι ήταν μόνο αυτό :lool: