# DSP lab9 part2

Apr 22nd, 2021
509
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
1. function [V] = aeq(inputPath, outputPath, Q, gain1, gain2, gain3, gain4, gain5)
2.
4.
5. Z = Y;
6.
7. Y_fft = fft(Y);
8. freqInitPlot = ((0:length(Y_fft) - 1) * FS) / length(Y_fft);
9.
10. figure;
11. plot(freqInitPlot, abs(Y_fft));
12. xlabel('Frequency (Hz)');
13. ylabel('Magnitude');
14. title('Initial FFT magnitude');
15.
16.
17. freqC1 = 60; % Hz
18. freqC2 = 230;
19. freqC3 = 910;
20. freqC4 = 4000;
21. freqC5 = 14000;
22.
23. freqL1 = freqC1 * ( (-1) / (2 * Q) + sqrt(1 + (1 / 4 * Q ^ 2)));
24. freqL2 = freqC2 * ( (-1) / (2 * Q) + sqrt(1 + (1 / 4 * Q ^ 2)));
25. freqL3 = freqC3 * ( (-1) / (2 * Q) + sqrt(1 + (1 / 4 * Q ^ 2)));
26. freqL4 = freqC4 * ( (-1) / (2 * Q) + sqrt(1 + (1 / 4 * Q ^ 2)));
27. freqL5 = freqC5 * ( (-1) / (2 * Q) + sqrt(1 + (1 / 4 * Q ^ 2)));
28.
29. freqU1 = freqC1 * ( 1 / (2 * Q) + sqrt(1 + (1 / 4 * Q ^ 2)));
30. freqU2 = freqC2 * ( 1 / (2 * Q) + sqrt(1 + (1 / 4 * Q ^ 2)));
31. freqU3 = freqC3 * ( 1 / (2 * Q) + sqrt(1 + (1 / 4 * Q ^ 2)));
32. freqU4 = freqC4 * ( 1 / (2 * Q) + sqrt(1 + (1 / 4 * Q ^ 2)));
33. freqU5 = freqC5 * ( 1 / (2 * Q) + sqrt(1 + (1 / 4 * Q ^ 2)));
34.
35. normFL1 = freqL1 / (FS / 2);
36. normFL2 = freqL2 / (FS / 2);
37. normFL3 = freqL3 / (FS / 2);
38. normFL4 = freqL4 / (FS / 2);
39. normFL5 = freqL5 / (FS / 2);
40.
41. normFU1 = freqU1 / (FS / 2);
42. normFU2 = freqU2 / (FS / 2);
43. normFU3 = freqU3 / (FS / 2);
44. normFU4 = freqU4 / (FS / 2);
45. normFU5 = freqU5 / (FS / 2);
46.
47. bw1 = freqC1 / Q;
48. bw2 = freqC2 / Q;
49. bw3 = freqC3 / Q;
50. bw4 = freqC4 / Q;
51. bw5 = freqC5 / Q;
52.
53. nOrder = 40;
54.
55. [filt1Num, filt1Den] = fir1(nOrder, normFU1, 'low');
56. [filt2Num, filt2Den] = fir1(nOrder, [normFL2 normFU2], 'bandpass');
57. [filt3Num, filt3Den] = fir1(nOrder, [normFL3 normFU3], 'bandpass');
58. [filt4Num, filt4Den] = fir1(nOrder, [normFL4 normFU4], 'bandpass');
59. [filt5Num, filt5Den] = fir1(nOrder, normFL5, 'high');
60.
61. filtered1 = filter(filt1Num, filt1Den, Y);
62. filtered2 = filter(filt2Num, filt2Den, Y);
63. filtered3 = filter(filt3Num, filt3Den, Y);
64. filtered4 = filter(filt4Num, filt4Den, Y);
65. filtered5 = filter(filt5Num, filt5Den, Y);
66.
67. filtG1 = filtered1 * 10 ^ (gain1 / 20);
68. filtG2 = filtered2 * 10 ^ (gain2 / 20);
69. filtG3 = filtered3 * 10 ^ (gain3 / 20);
70. filtG4 = filtered4 * 10 ^ (gain4 / 20);
71. filtG5 = filtered5 * 10 ^ (gain5 / 20);
72.
73. Z = filtG1 + filtG2 + filtG3 + filtG4 + filtG5;
74.
75. Z_fft = fft(Z);
76. freqInitPlot = ((0:length(Z_fft) - 1) * FS) / length(Z_fft);
77.
78. figure;
79. plot(freqInitPlot, abs(Z_fft));
80. xlabel('Frequency (Hz)');
81. ylabel('Magnitude');
82. title('Output FFT magnitude');
83.
84. audiowrite(outputPath, Z, FS);
85.
86. end
RAW Paste Data