- 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

4 views (last 30 days)

Show older comments

##### 0 Comments

### Answers (1)

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 :

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

Mathieu NOE
on 2 Feb 2022

hello

do you mind accepting my answer (if it fullfills your expectations) ?

tx

### See Also

### Community Treasure Hunt

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

Start Hunting!