hello
as an appetizer , please read the attached txt file and see how it's implemented in the m file (filters.m)
for your specific request , the code below will do these tasks :
- import the audio file signal
- design the filter
- apply the filter to the audio signal
- as a bonus , I added a last section to show how a fft can be used on signals (the raw and the filtered) to compute the transfer function of the filter applied on the signals without having any a priori information about the filter ; you can see that the bode plot generated in the filter design phase is the same as the one obtained by fft on the audio signals
hope this helps
clc
clearvars
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% load signal
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[signal,Fs] = audioread('noisy2_group2.wav');
[samples,channels] = size(signal);
signal = signal(:,1); % select first channel (example) if the audio file is multichannel
% create time vector
dt = 1/Fs;
time = (0:samples-1)*dt;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% highShelf: H(s) = A * (A*s^2 + (sqrt(A)/Q)*s + 1)/(s^2 + (sqrt(A)/Q)*s + A)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
A = 4;
Q = 3;
f0 = 100;
w0 = 2*pi*f0/Fs
alpha = sin(w0)/(2*Q);
b0 = A*( (A+1) + (A-1)*cos(w0) + 2*sqrt(A)*alpha );
b1 = -2*A*( (A-1) + (A+1)*cos(w0) );
b2 = A*( (A+1) + (A-1)*cos(w0) - 2*sqrt(A)*alpha );
a0 = (A+1) - (A-1)*cos(w0) + 2*sqrt(A)*alpha;
a1 = 2*( (A-1) - (A+1)*cos(w0) );
a2 = (A+1) - (A-1)*cos(w0) - 2*sqrt(A)*alpha;
B = [b0 b1 b2];
A = [a0 a1 a2];
% Bode plot
freq = logspace(1,(log10(Fs/2.5)),300);
[g1,p1] = dbode(B,A,1/Fs,2*pi*freq);
g1db = 20*log10(g1);
figure(1);
subplot(2,1,1),semilogx(freq,g1db);
title('highShelf: H(s) = A * (A*s^2 + (sqrt(A)/Q)*s + 1)/(s^2 + (sqrt(A)/Q)*s + A)');
ylabel('dB ');
subplot(2,1,2),semilogx(freq,p1);
xlabel('Fréquence (Hz)');
ylabel(' phase (°)')
% apply filter to signal
signal_filtered = filter(B,A,signal);
% if you want to export to wav format , you have to normalize the
% amplitude of the filtered signal into +/- 1 range
mm = max(abs(signal_filtered));
if mm>=1
signal_filtered = signal_filtered./mm;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display 1 : time domain plot
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(2),
subplot(2,1,1),plot(time,signal,'b');
title(['Time plot / Fs = ' num2str(Fs) ' Hz ']);
ylabel('Amplitude');
legend('original');
subplot(2,1,2),plot(time,signal_filtered,'r');
xlabel('Time (s)');ylabel('Amplitude');
legend('filtered');
nfft = 512*8;
window = hanning(nfft);
noverlap = round(0.75*nfft);
[Txy,Cxy,Fv] = tfe_and_coh(signal,signal_filtered,nfft,Fs,window,noverlap);
figure(3)
subplot(3,1,1)
semilogx(Fv, 20*log10(abs(Txy)))
title('Amplitude')
ylabel('dB')
subplot(3,1,2)
semilogx(Fv, angle(Txy)*180/pi)
title('Phase')
ylabel('°')
subplot(3,1,3)
semilogx(Fv, Cxy)
title('Coherence')
%%%%%%%%%%%%%%%%%%%%%%%
function [Txy,Cxy,f] = tfe_and_coh(x,y,nfft,Fs,window,noverlap)
% Transfer Function and Coherence Estimate
% compute PSD and CSD
window = window(:);
n = length(x); % Number of data points
nwind = length(window); % length of window
if n < nwind % zero-pad x , y if length is less than the window length
x(nwind)=0;
y(nwind)=0;
n=nwind;
end
x = x(:); % Make sure x is a column vector
y = y(:); % Make sure y is a column vector
k = fix((n-noverlap)/(nwind-noverlap)); % Number of windows
% (k = fix(n/nwind) for noverlap=0)
index = 1:nwind;
Pxx = zeros(nfft,1);
Pyy = zeros(nfft,1);
Pxy = zeros(nfft,1);
for i=1:k
xw = window.*x(index);
yw = window.*y(index);
index = index + (nwind - noverlap);
Xx = fft(xw,nfft);
Yy = fft(yw,nfft);
Xx2 = abs(Xx).^2;
Yy2 = abs(Yy).^2;
Xy2 = Yy.*conj(Xx);
Pxx = Pxx + Xx2;
Pyy = Pyy + Yy2;
Pxy = Pxy + Xy2;
end
% Select first half
if ~any(any(imag([x y])~=0)), % if x and y are not complex
if rem(nfft,2), % nfft odd
select = [1:(nfft+1)/2];
else
select = [1:nfft/2+1]; % include DC AND Nyquist
end
Pxx = Pxx(select);
Pyy = Pyy(select);
Pxy = Pxy(select);
else
select = 1:nfft;
end
Txy = Pxy ./ Pxx; % transfer function estimate
Cxy = (abs(Pxy).^2)./(Pxx.*Pyy); % coherence function estimate
f = (select - 1)'*Fs/nfft;
end
