tsounakis

2.4 Signal Processing

Nov 28th, 2021 (edited)
332
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. %Αθανάσιος Τσουνάκης 2021
  2. %Εκφώνηση 2.4:
  3. %Να σχεδιαστούν δύο φίλτρα δεύτερης και τρίτης τάξης ψηφιακό Butterworth
  4. %ζωνοπερατό φίλτρο με τις εξής προδιαγραφές:
  5. % >συχνότητες αποκοπής 80 Hz και 150 Hz
  6. % >συχνότητα δειγματοληψίας 48 KHz
  7. %Στη συνέχεια να σχεδιαστούν στο ίδιο γράφημα οι καμπύλες πλάτους και φάσης
  8. %των δύο φίλτρων διαφορετικής τάξης.
  9.  
  10. function [] = f()                                      
  11.     clc;
  12.     clear;
  13.     syms s
  14.        
  15.     Fs = 48*10^3;
  16.    
  17.     % βήμα 1 και 2: εύρεση επιθυμητών συχνοτήτων και (ψηφιακών) και
  18.     % αντιστάθμιση της παραμόρφωσης που θα προέλθει λόγω του διγραμμικού
  19.     % μετασχηματισμού (στην συνάρτηση digitalAngular)
  20.     [cutoffL, cutoffH] = digitalAngular(Fs, 0.08*10^3, 0.15*10^3);
  21.    
  22.     % βήμα 3α: επιλογή των πρωτότυπων συναρτήσεων μεταφοράς πρώτης και
  23.     % δεύτερης τάξης Butterworth φίλτρου
  24.     transferButterworthFirst(s) = 1/(s+1);
  25.     transferButterworthSec(s) = 1/(s^2+1.4142*s+1);
  26.    
  27.     % υπόλοιπα βήματα 3β+
  28.     plotTF(cutoffH, cutoffL, transferButterworthFirst, transferButterworthSec, Fs);
  29.    
  30. end
  31.  
  32. function plotTF(cutoffH, cutoffL, trff, trfs, Fs)
  33.     syms s
  34.     syms z
  35.    
  36.     % ορισμός της κεντρικής συχνότητας resonant frequency (ω_c)
  37.     resonantFreq = (cutoffH*cutoffL)^(0.5);
  38.    
  39.     % κατασκευή της συνάρτησης μεταφοράς των φίλτρων (δεύτερης και τρίτης
  40.     % τάξης ψηφιακών Butterworth ζωνοπερατών φίλτρων)
  41.     transferBPF(s) = simplify(trff((s^2+resonantFreq^2)/(s*(cutoffH-cutoffL))))
  42.     transferBPS(s) = simplify(trfs((s^2+resonantFreq^2)/(s*(cutoffH-cutoffL))))
  43.        
  44.     % μετασχηματισμός των δύο συναρτήσεων από συμβολικές εκφράσεις ως προς
  45.     % s σε αριθμητικές εκφράσεις (πρακτικό ζήτημα Matlab)
  46.     transferSPlaneF = syms2numeric(transferBPF(s));
  47.     [transferSPlaneDenF, transferSPlaneNumF] = tf2string(transferSPlaneF);
  48.     transferSPlaneS = syms2numeric(transferBPS(s));
  49.     [transferSPlaneDenS, transferSPlaneNumS] = tf2string(transferSPlaneS);
  50.    
  51.     tSPlaneStr = sprintf('$%s=\\displaystyle\\frac{%s}{%s}$ ', 'H(\omega)',...
  52.         transferSPlaneNumF, transferSPlaneDenF);
  53.    
  54.     % εμφάνιση των γραφημάτων (κρουστική απόκριση και συχνοτική απόκριση)
  55.     figure(1)
  56.     subplot(2,1,1);
  57.     bodemag(transferSPlaneF, transferSPlaneS); title(sprintf('\\textbf{Magnitude Diagram, ${\\omega_c = %s}$ ${rad/s }$ }'...
  58.         , num2str(resonantFreq,'%.3f')), 'Interpreter', 'latex'); xlabel('$Frequency$', 'Interpreter', 'latex');
  59.     ylabel('$Phase$', 'Interpreter', 'latex'); xlabel('$Frequency$', 'Interpreter', 'latex');
  60.     ylabel('$Magnitude$', 'Interpreter', 'latex');
  61.     legend('Second-Order Filter', 'Third-Order Filter', '', 'Resonant Frequency','High cut-off Frequency');
  62.     xline(cutoffL, 'r-.', 'DisplayName', 'Lower Cut-Off Frequency');
  63.     xline(resonantFreq, 'b--', 'DisplayName', 'Resonant Frequency');
  64.     xline(cutoffH, 'r-.', 'DisplayName', 'Higher Cut-Off Frequency');
  65.     ylim([-140 20]);
  66.     grid on;
  67.    
  68.     figure(1)
  69.     subplot(2,1,2);
  70.     h = bodeplot(transferSPlaneF, transferSPlaneS);
  71.     setoptions(h,'MagVisible','off'); t = title(sprintf('\\textbf{Phase Diagram, ${\\omega_c = %s}$ ${rad/s }$ }', num2str(resonantFreq,'%.3f')), 'Interpreter', 'latex'); xlabel('$Frequency$', 'Interpreter', 'latex');
  72.     ylabel('$Phase$', 'Interpreter', 'latex');
  73.     legend('Second-Order DBBF','Third-Order DBBF', 'Low cut-off Frequency', 'Resonant Frequency','High cut-off Frequency');
  74.     xline(cutoffL, 'r-.', 'DisplayName', 'Lower Cut-Off Frequency');
  75.     xline(resonantFreq, 'b--', 'DisplayName', 'Resonant Frequency');
  76.     xline(cutoffH, 'r-.', 'DisplayName', 'Higher Cut-Off Frequency');
  77.     grid on;
  78.    
  79.     % υπολογισμός των συναρτήσεων μεταφοράς των φίλτρων στο πεδίο-z
  80.     transferZPlaneF = transferBPF(2*((z-1)/(z+1))*(Fs))
  81.     transferZPlaneS = transferBPS(2*((z-1)/(z+1))*(Fs))
  82.    
  83. end
  84.  
  85. function [numericTransfer] = syms2numeric(G)
  86.     [symNum,symDen] = numden(G);
  87.     TFnum = sym2poly(symNum);    
  88.     TFden = sym2poly(symDen);    
  89.     numericTransfer = tf(TFnum,TFden);
  90. end
  91.  
  92. function [cutoffL, cutoffH] = digitalAngular(Fs, cutoffLA, cutoffLH)
  93.     T = 1 / Fs;
  94.     cutoffL = 2*tan(2*pi*cutoffLA*T/2)/T;
  95.     cutoffH = 2*tan(2*pi*cutoffLH*T/2)/T;
  96. end
  97.  
  98. function [char_den, char_num] = tf2string(input_tf) %modified version of the
  99.     syms s                                          %function made by vittorio88 (stackoverflow)
  100.        
  101.     sym_num=poly2sym(input_tf.num{:},s);
  102.     sym_num=vpa(sym_num, 4);
  103.     char_num=char(sym_num);
  104.  
  105.     sym_den=poly2sym(input_tf.den{:},s);
  106.     sym_den=vpa(sym_den, 4);
  107.     char_den=char(sym_den);
  108.  
  109.     s=tf('s');
  110.  
  111. end
  112.  
RAW Paste Data