DSP lab9 part2

function [V] = aeq(inputPath, outputPath, Q, gain1, gain2, gain3, gain4, gain5)

[Y, FS] = audioread(inputPath);

Z = Y;

Y_fft = fft(Y);
freqInitPlot = ((0:length(Y_fft) - 1) * FS) / length(Y_fft);

figure;
plot(freqInitPlot, abs(Y_fft));
xlabel('Frequency (Hz)');
ylabel('Magnitude');
title('Initial FFT magnitude');


freqC1 = 60; % Hz
freqC2 = 230;
freqC3 = 910;
freqC4 = 4000;
freqC5 = 14000;

freqL1 = freqC1 * ( (-1) / (2 * Q) + sqrt(1 + (1 / 4 * Q ^ 2)));
freqL2 = freqC2 * ( (-1) / (2 * Q) + sqrt(1 + (1 / 4 * Q ^ 2)));
freqL3 = freqC3 * ( (-1) / (2 * Q) + sqrt(1 + (1 / 4 * Q ^ 2)));
freqL4 = freqC4 * ( (-1) / (2 * Q) + sqrt(1 + (1 / 4 * Q ^ 2)));
freqL5 = freqC5 * ( (-1) / (2 * Q) + sqrt(1 + (1 / 4 * Q ^ 2)));

freqU1 = freqC1 * ( 1 / (2 * Q) + sqrt(1 + (1 / 4 * Q ^ 2)));
freqU2 = freqC2 * ( 1 / (2 * Q) + sqrt(1 + (1 / 4 * Q ^ 2)));
freqU3 = freqC3 * ( 1 / (2 * Q) + sqrt(1 + (1 / 4 * Q ^ 2)));
freqU4 = freqC4 * ( 1 / (2 * Q) + sqrt(1 + (1 / 4 * Q ^ 2)));
freqU5 = freqC5 * ( 1 / (2 * Q) + sqrt(1 + (1 / 4 * Q ^ 2)));

normFL1 = freqL1 / (FS / 2);
normFL2 = freqL2 / (FS / 2);
normFL3 = freqL3 / (FS / 2);
normFL4 = freqL4 / (FS / 2);
normFL5 = freqL5 / (FS / 2);

normFU1 = freqU1 / (FS / 2);
normFU2 = freqU2 / (FS / 2);
normFU3 = freqU3 / (FS / 2);
normFU4 = freqU4 / (FS / 2);
normFU5 = freqU5 / (FS / 2);

bw1 = freqC1 / Q;
bw2 = freqC2 / Q;
bw3 = freqC3 / Q;
bw4 = freqC4 / Q;
bw5 = freqC5 / Q;

nOrder = 40;

[filt1Num, filt1Den] = fir1(nOrder, normFU1, 'low');
[filt2Num, filt2Den] = fir1(nOrder, [normFL2 normFU2], 'bandpass');
[filt3Num, filt3Den] = fir1(nOrder, [normFL3 normFU3], 'bandpass');
[filt4Num, filt4Den] = fir1(nOrder, [normFL4 normFU4], 'bandpass');
[filt5Num, filt5Den] = fir1(nOrder, normFL5, 'high');

filtered1 = filter(filt1Num, filt1Den, Y);
filtered2 = filter(filt2Num, filt2Den, Y);
filtered3 = filter(filt3Num, filt3Den, Y);
filtered4 = filter(filt4Num, filt4Den, Y);
filtered5 = filter(filt5Num, filt5Den, Y);

filtG1 = filtered1 * 10 ^ (gain1 / 20);
filtG2 = filtered2 * 10 ^ (gain2 / 20);
filtG3 = filtered3 * 10 ^ (gain3 / 20);
filtG4 = filtered4 * 10 ^ (gain4 / 20);
filtG5 = filtered5 * 10 ^ (gain5 / 20);

Z = filtG1 + filtG2 + filtG3 + filtG4 + filtG5;

Z_fft = fft(Z);
freqInitPlot = ((0:length(Z_fft) - 1) * FS) / length(Z_fft);

figure;
plot(freqInitPlot, abs(Z_fft));
xlabel('Frequency (Hz)');
ylabel('Magnitude');
title('Output FFT magnitude');

audiowrite(outputPath, Z, FS);

end