i am trying to code a high shelf filter with an input(.wav) audio signal . I am a student im doing this for my final year project....kindly professionals give me solution

4 views (last 30 days)
Khader Mydeen
Khader Mydeen on 16 Jan 2022
Commented: Mathieu NOE on 2 Feb 2022
i am trying to code a high shelf filter with an input(.wav) audio signal . I am a student im doing this for my final year project....kindly professionals give me solution...........i dont even know basics....Kindly help

Answers (1)

Mathieu NOE
Mathieu NOE on 17 Jan 2022
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
  4 Comments

Sign in to comment.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!