- 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
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
3 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
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
0 comentarios
Respuestas (1)
Mathieu NOE
el 17 de En. de 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 :
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 comentarios
Mathieu NOE
el 2 de Feb. de 2022
hello
do you mind accepting my answer (if it fullfills your expectations) ?
tx
Ver también
Categorías
Más información sobre Audio Processing Algorithm Design en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!