trouble making a notchfilter to remove a tone out of a wavefile

Question

Jacob 2023-11-29

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2053427-trouble-making-a-notchfilter-to-remove-a-tone-out-of-a-wavefile

评论： Mathieu NOE 2023-12-1

I am trying to write a matlab code to design a notchfilter to filter a tone out of a wav.file but I am having trouble getting the filter to work this code below is supposed to determine the frequencies of the audio in the wav.file and filter the tone out of the sunshinesquare.wav file. however when I go to run the program it plays the audio file with the tone in it still. I can't figure out what is wrong with it any suggestions/help is greatly appreciated.

My code is below

[sound, fs] = audioread('SunshineSquare.wav');
soundsc(sound,fs)
figure(1)
specgram(sound,fs)
p1 = 2 * pi; 
bb = poly([exp(1*p1/4) exp(-1*p1/4) exp(1*p1/2) exp(-1*p1/2)]);
aa = poly(0.9 * [exp(1*p1/4) exp(-1*p1/4) exp(1*p1/2) exp(-1*p1/2)]);
figure;
zplane(bb, aa);
title('Poles and Zeros Plot');
Ww = -p1:p1/100:p1;
HH1 = freqz(bb, 1, Ww);
HH2 = freqz(bb, aa, Ww);
figure;
plot(Ww, abs(HH1), 'b', 'LineWidth', 2);
hold on;
plot(Ww, abs(HH2), 'r', 'LineWidth', 2);
hold off;
legend('Zeros Only', 'Poles and Zeros');
xlabel('Frequency');
ylabel('Magnitude');
title('Frequency Response');
audiowrite('processed_SunshineSquare.wav', y_processed, Fs_processed);
soundsc(y, Fs);

2 个评论
显示无隐藏无

Walter Roberson 2023-11-29

cross-reference https://www.mathworks.com/matlabcentral/answers/2053422-trouble-making-a-firfilter-to-remove-a-tone-out-of-a-wavefile?s_tid=srchtitle

Walter Roberson 2023-11-29

在 MATLAB Online 中打开

audiowrite('processed_SunshineSquare.wav', y_processed, Fs_processed);
soundsc(y, Fs);

None of those variables is defined in your code. Not even Fs.

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Mathieu NOE 2023-11-29

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2053427-trouble-making-a-notchfilter-to-remove-a-tone-out-of-a-wavefile#answer_1362237

编辑：Mathieu NOE 2023-11-29

在 MATLAB Online 中打开

hello

adapt this code to your needs

you can apply multiples notches , simply create the list of frequencies here

fc_notch = [9888; 10896]; % notch center frequencies

cheers

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% select filters to apply
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
option_notch = 1; % 0 = without notch filter , 1 = with notch filter
    fc_notch = [9888; 10896];       % notch center frequencies
    p = 0.98;  % bandwith parameter (closer to 1 reduces the BW)
option_BPF = 1; % 0 = without bandpass filter , 1 = with bandpass filter
    fc = 1e4;      % Hz
    bandwith = 1e4;      % Hz
    N_bpf = 4;           % filter order
     
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% options 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% if you are dealing with acoustics, you may wish to have A weighted
% spectrums 
% option_w = 0 : linear spectrum (no weighting dB (L) )
% option_w = 1 : A weighted spectrum (dB (A) )
option_w = 0;     
% spectrogram dB scale
spectrogram_dB_scale = 80;  % the lowest displayed level is XX dB below the max level
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% load signal
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[signal,Fs] = audioread('original.wav');
dt = 1/Fs;
[samples,channels] = size(signal);
% select channel (if needed)
channels = 1;
signal = signal(:,channels);
% time vector 
time = (0:samples-1)*dt;
%% decimate (if needed)
% NB : decim = 1 will do nothing (output = input)
decim = 0;
if decim>1
    Fs = Fs/decim;
    for ck = 1:channels
    newsignal(:,ck) = decimate(signal(:,ck),decim);
    end
   signal = newsignal;
end
signal_filtered = signal;
samples = length(signal);
time = (0:samples-1)*1/Fs;
%% notch filter section %%%%%%
% y(n)=G*[x(n)-2*cos(w0)*x(n-1)+x(n-2)]+[2*p cos(w0)*y(n-1)-p^2 y(n-2)]
if option_notch ~= 0
    for ck =1:numel(fc_notch)
        w0 = 2*pi*fc_notch(ck)/Fs;
        % digital notch (IIR)
        num1z=[1 -2*cos(w0) 1];
        den1z=[1 -2*p*cos(w0) p^2];
        % now let's filter the signal 
        signal_filtered = filter(num1z,den1z,signal_filtered);
    end
end
%% band pass  filter section %%%%%%
if option_BPF ~= 0
    f_low = fc-0.5*bandwith;
    f_high = fc+0.5*bandwith;
    [b,a] = butter(N_bpf,2/Fs*[f_low f_high]);
    signal_filtered = filter(b,a,signal_filtered);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FFT parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
NFFT = 1000;    % 
OVERLAP = 0.75;  % percentage of overlap
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display 1 : time domain plot
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(1),
subplot(211),plot(time,signal,'b');grid on
title(['Time plot  / Fs = ' num2str(Fs) ' Hz / raw data ']);
xlabel('Time (s)');ylabel('Amplitude');
subplot(212),plot(time,signal_filtered,'r');grid on
title(['Time plot  / Fs = ' num2str(Fs) ' Hz / filtered data ']);
xlabel('Time (s)');ylabel('Amplitude');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display 2 : averaged FFT spectrum
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[freq, spectrum_raw] = myfft_peak(signal,Fs,NFFT,OVERLAP);
[freq, spectrum_filtered] = myfft_peak(signal_filtered,Fs,NFFT,OVERLAP);
figure(2),plot(freq,20*log10(spectrum_raw),'b',freq,20*log10(spectrum_filtered),'r');grid on
df = freq(2)-freq(1); % frequency resolution 
title(['Averaged FFT Spectrum  / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(df,3) ' Hz ']);
xlabel('Frequency (Hz)');ylabel('Amplitude (dB (L))');
legend('raw','filtered');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display 3 : time / frequency analysis : spectrogram demo
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for ck = 1:channels
    [sg,fsg,tsg] = specgram(signal(:,ck),NFFT,Fs,hanning(NFFT),floor(NFFT*OVERLAP));  
    % FFT normalisation and conversion amplitude from linear to dB (peak)
    sg_dBpeak = 20*log10(abs(sg))+20*log10(2/length(fsg));     % NB : X=fft(x.*hanning(N))*4/N; % hanning only
    % apply A weigthing if needed
    if option_w == 1
        pondA_dB = pondA_function(fsg);
        sg_dBpeak = sg_dBpeak+(pondA_dB*ones(1,size(sg_dBpeak,2)));
        my_title = ('Spectrogram (dB (A))');
    else
        my_title = ('Spectrogram (dB (L))');
    end
    % saturation of the dB range :  the lowest displayed level is XX dB below the max level
    min_disp_dB = round(max(max(sg_dBpeak))) - spectrogram_dB_scale;
    sg_dBpeak(sg_dBpeak<min_disp_dB) = min_disp_dB;
    % plots spectrogram
    figure(2+ck);
    imagesc(tsg,fsg,sg_dBpeak);colormap('jet');
    axis('xy');colorbar('vert');grid on
    df = fsg(2)-fsg(1); % freq resolution 
    title([my_title ' / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(df,3) ' Hz / Channel : ' num2str(ck)]);
    xlabel('Time (s)');ylabel('Frequency (Hz)');
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function pondA_dB = pondA_function(f)
	% dB (A) weighting curve
	n = ((12200^2*f.^4)./((f.^2+20.6^2).*(f.^2+12200^2).*sqrt(f.^2+107.7^2).*sqrt(f.^2+737.9^2)));
	r = ((12200^2*1000.^4)./((1000.^2+20.6^2).*(1000.^2+12200^2).*sqrt(1000.^2+107.7^2).*sqrt(1000.^2+737.9^2))) * ones(size(f));
	pondA = n./r;
	pondA_dB = 20*log10(pondA(:));
end
function  [freq_vector,fft_spectrum] = myfft_peak(signal, Fs, nfft, Overlap)
% FFT peak spectrum of signal  (example sinus amplitude 1   = 0 dB after fft).
% Linear averaging
%   signal - input signal, 
%   Fs - Sampling frequency (Hz).
%   nfft - FFT window size
%   Overlap - buffer percentage of overlap % (between 0 and 0.95)
[samples,channels] = size(signal);
% fill signal with zeros if its length is lower than nfft
if samples<nfft
    s_tmp = zeros(nfft,channels);
    s_tmp((1:samples),:) = signal;
    signal = s_tmp;
    samples = nfft;
end
% window : hanning
window = hanning(nfft);
window = window(:);
%    compute fft with overlap 
 offset = fix((1-Overlap)*nfft);
 spectnum = 1+ fix((samples-nfft)/offset); % Number of windows
%     % for info is equivalent to : 
%     noverlap = Overlap*nfft;
%     spectnum = fix((samples-noverlap)/(nfft-noverlap));	% Number of windows
    % main loop
    fft_spectrum = 0;
    for i=1:spectnum
        start = (i-1)*offset;
        sw = signal((1+start):(start+nfft),:).*(window*ones(1,channels));
        fft_spectrum = fft_spectrum + (abs(fft(sw))*4/nfft);     % X=fft(x.*hanning(N))*4/N; % hanning only 
    end
    fft_spectrum = fft_spectrum/spectnum; % to do linear averaging scaling
% one sidded fft spectrum  % Select first half 
    if rem(nfft,2)    % nfft odd
        select = (1:(nfft+1)/2)';
    else
        select = (1:nfft/2+1)';
    end
fft_spectrum = fft_spectrum(select,:);
freq_vector = (select - 1)*Fs/nfft;
end

13 个评论
显示 11更早的评论隐藏 11更早的评论

Mathieu NOE 2023-11-30

编辑：Mathieu NOE 2023-11-30

在 MATLAB Online 中打开

@Jacob

try this

the notch filters are now tuned to match the 3 tones (the last one is at 4725 Hz)

no band pass filter is activated here (option_BPF = 0;)

always double check that filters must have their corner / notch frequencies below Fs/2 (aka Nyquist freq) otherwise you get the above error message

also I have modified a bit my code so that each channel of data gets one figure with raw signal spectrogram on top and filtered signal below

NB that the supplied file is mono and not 2 channels but that is not an issue in itself , just fyi

all the best

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% select audio file 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
audiofile = 'SunshineSquare.wav';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% select filters to apply
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
option_notch = 1; % 0 = without notch filter , 1 = with notch filter
    fc_notch = 4725*[1/3 2/3 1];       % notch center frequencies
    p = 0.95;  % bandwith parameter (closer to 1 reduces the BW)
option_BPF = 0; % 0 = without bandpass filter , 1 = with bandpass filter
    fc = 1e4;      % Hz
    bandwith = 1e4;      % Hz
    N_bpf = 4;           % filter order
     
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% options 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% if you are dealing with acoustics, you may wish to have A weighted
% spectrums 
% option_w = 0 : linear spectrum (no weighting dB (L) )
% option_w = 1 : A weighted spectrum (dB (A) )
option_w = 0;     
% spectrogram dB scale
spectrogram_dB_scale = 80;  % the lowest displayed level is XX dB below the max level
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% load signal
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[signal,Fs] = audioread(audiofile);
[samples,channels] = size(signal);
dt = 1/Fs;
% % select channel (if needed)
% channels = 1:2;
% signal = signal(:,channels);
% [samples,channels] = size(signal);
% time vector 
time = (0:samples-1)*dt;
%% decimate (if needed)
% NB : decim = 1 will do nothing (output = input)
decim = 0;
if decim>1
    Fs = Fs/decim;
    for ck = 1:channels
    newsignal(:,ck) = decimate(signal(:,ck),decim);
    end
   signal = newsignal;
end
signal_filtered = signal;
samples = length(signal);
time = (0:samples-1)*1/Fs;
%% notch filter section %%%%%%
% y(n)=G*[x(n)-2*cos(w0)*x(n-1)+x(n-2)]+[2*p cos(w0)*y(n-1)-p^2 y(n-2)]
if option_notch ~= 0
    for ck =1:numel(fc_notch)
        w0 = 2*pi*fc_notch(ck)/Fs;
        % digital notch (IIR)
        num1z=[1 -2*cos(w0) 1];
        den1z=[1 -2*p*cos(w0) p^2];
        % now let's filter the signal 
        signal_filtered = filter(num1z,den1z,signal_filtered);
    end
end
%% band pass  filter section %%%%%%
if option_BPF ~= 0
    f_low = fc-0.5*bandwith;
    f_high = fc+0.5*bandwith;
    [b,a] = butter(N_bpf,2/Fs*[f_low f_high]);
    signal_filtered = filter(b,a,signal_filtered);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FFT parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
NFFT = 1000;    % 
OVERLAP = 0.75;  % percentage of overlap
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display 1 : time domain plot
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for ck = 1:channels
    figure(ck),
    subplot(211),plot(time,signal(:,ck),'b');grid on
    title(['Time plot  / Fs = ' num2str(Fs) ' Hz / raw data / Channel : ' num2str(ck)]);
    xlabel('Time (s)');ylabel('Amplitude');
    subplot(212),plot(time,signal_filtered(:,ck),'r');grid on
    title(['Time plot  / Fs = ' num2str(Fs) ' Hz / filtered data / Channel : ' num2str(ck)]);
    xlabel('Time (s)');ylabel('Amplitude');
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display 2 : averaged FFT spectrum
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[freq, spectrum_raw] = myfft_peak(signal,Fs,NFFT,OVERLAP);
[freq, spectrum_filtered] = myfft_peak(signal_filtered,Fs,NFFT,OVERLAP);
df = freq(2)-freq(1); % frequency resolution 
for ck = 1:channels
    figure(channels+ck),
    plot(freq,20*log10(spectrum_raw(:,ck)),'b',freq,20*log10(spectrum_filtered(:,ck)),'r');grid on
    title(['Averaged FFT Spectrum  / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(df,3) ' Hz / Channel : ' num2str(ck)]);
    xlabel('Frequency (Hz)');ylabel('Amplitude (dB (L))');
    legend('raw','filtered');
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display 3 : time / frequency analysis : spectrogram demo
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for ck = 1:channels
    [sg,fsg,tsg] = specgram(signal(:,ck),NFFT,Fs,hanning(NFFT),floor(NFFT*OVERLAP));  
    [sgf,fsgf,tsgf] = specgram(signal_filtered(:,ck),NFFT,Fs,hanning(NFFT),floor(NFFT*OVERLAP));  
    % FFT normalisation and conversion amplitude from linear to dB (peak)
    sg_dBpeak = 20*log10(abs(sg))+20*log10(2/length(fsg));     % NB : X=fft(x.*hanning(N))*4/N; % hanning only
    sgf_dBpeak = 20*log10(abs(sgf))+20*log10(2/length(fsg));     % NB : X=fft(x.*hanning(N))*4/N; % hanning only
    % apply A weigthing if needed
    if option_w == 1
        pondA_dB = pondA_function(fsg);
        sg_dBpeak = sg_dBpeak+(pondA_dB*ones(1,size(sg_dBpeak,2)));
        sgf_dBpeak = sgf_dBpeak+(pondA_dB*ones(1,size(sgf_dBpeak,2)));
        my_title = ('Spectrogram (dB (A))');
    else
        my_title = ('Spectrogram (dB (L))');
    end
    % saturation of the dB range :  the lowest displayed level is XX dB below the max level
    min_disp_dB = round(max(max(sg_dBpeak))) - spectrogram_dB_scale;
    sg_dBpeak(sg_dBpeak<min_disp_dB) = min_disp_dB;
    sgf_dBpeak(sgf_dBpeak<min_disp_dB) = min_disp_dB;
    % plots spectrogram
    figure(2*channels+ck);
    subplot(2,1,1),imagesc(tsg,fsg,sg_dBpeak);colormap('jet');
    axis('xy');colorbar('vert');grid on
    df = fsg(2)-fsg(1); % freq resolution 
    title([my_title ' Raw Signal / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(df,3) ' Hz / Channel : ' num2str(ck)]);
    xlabel('Time (s)');ylabel('Frequency (Hz)');
    subplot(2,1,2),imagesc(tsgf,fsgf,sgf_dBpeak);colormap('jet');
    axis('xy');colorbar('vert');grid on
    title([my_title ' / Filtered Signal  / Channel : ' num2str(ck)]);
    xlabel('Time (s)');ylabel('Frequency (Hz)');
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function pondA_dB = pondA_function(f)
	% dB (A) weighting curve
	n = ((12200^2*f.^4)./((f.^2+20.6^2).*(f.^2+12200^2).*sqrt(f.^2+107.7^2).*sqrt(f.^2+737.9^2)));
	r = ((12200^2*1000.^4)./((1000.^2+20.6^2).*(1000.^2+12200^2).*sqrt(1000.^2+107.7^2).*sqrt(1000.^2+737.9^2))) * ones(size(f));
	pondA = n./r;
	pondA_dB = 20*log10(pondA(:));
end
function  [freq_vector,fft_spectrum] = myfft_peak(signal, Fs, nfft, Overlap)
% FFT peak spectrum of signal  (example sinus amplitude 1   = 0 dB after fft).
% Linear averaging
%   signal - input signal, 
%   Fs - Sampling frequency (Hz).
%   nfft - FFT window size
%   Overlap - buffer percentage of overlap % (between 0 and 0.95)
[samples,channels] = size(signal);
% fill signal with zeros if its length is lower than nfft
if samples<nfft
    s_tmp = zeros(nfft,channels);
    s_tmp((1:samples),:) = signal;
    signal = s_tmp;
    samples = nfft;
end
% window : hanning
window = hanning(nfft);
window = window(:);
%    compute fft with overlap 
 offset = fix((1-Overlap)*nfft);
 spectnum = 1+ fix((samples-nfft)/offset); % Number of windows
%     % for info is equivalent to : 
%     noverlap = Overlap*nfft;
%     spectnum = fix((samples-noverlap)/(nfft-noverlap));	% Number of windows
    % main loop
    fft_spectrum = 0;
    for i=1:spectnum
        start = (i-1)*offset;
        sw = signal((1+start):(start+nfft),:).*(window*ones(1,channels));
        fft_spectrum = fft_spectrum + (abs(fft(sw))*4/nfft);     % X=fft(x.*hanning(N))*4/N; % hanning only 
    end
    fft_spectrum = fft_spectrum/spectnum; % to do linear averaging scaling
% one sidded fft spectrum  % Select first half 
    if rem(nfft,2)    % nfft odd
        select = (1:(nfft+1)/2)';
    else
        select = (1:nfft/2+1)';
    end
fft_spectrum = fft_spectrum(select,:);
freq_vector = (select - 1)*Fs/nfft;
end

Mathieu NOE 2023-11-30

编辑：Mathieu NOE 2023-11-30

在 MATLAB Online 中打开

env_secant.m

hello again

from this point, this is just extra bonus , for my own fun

as you probably noticed while listening to the filtered sound, even though we succeeded to remove the 3 main disturbing frequencies (tx to the notch filter) , there are still some transients ("klicks") to be heared in the last portion of the signal

my new task that I assigned to myself (and accepted) was to figure out how to reduce these artifacts.

I put some code together , but as I progressed , I noticed that this new function would even be more efficient if the signal would not have those staircase effect (I suspect the person who created the noisy signal put not only extra tones but also positive and negative offsets - and we want to remove them)

that can be easily achieved if you add some highpass or bandpass filtering

good news, this is already available in the code already provided

simply change this portion of the code to this :

option_BPF = 1; % 0 = without bandpass filter , 1 = with bandpass filter
f_low = 20;
f_high = 5000;
N_bpf = 2;           % filter order

and now you have a signal that doesn't have those offsets anymore (the mean value remains always zero whatever the buffer length of signal you consider)

then we can start the last mile code :

(you can simply copy paste that to a new sheet if you don't want to mess up too much the main code)

you need env_secant.m which is in attachment

have fun !

% let's find out how we can make transients visible
ag = abs(gradient(signal_filtered)); % absolute of (gradient of signal) 
[env] = env_secant(time', ag, 120, 'top'); % secant envelope of ag (see function attached)
% figure(10),
% plot(time,ag,time,env)  % for your info
agm = sqrt(mean(ag.^2)); % rms value of ag => will be used to define a good threshold (line below)
TF = (env>12*agm); % find time indexes of  transients ("klicks")
% create copy signal and nullify it when the envelope of the abs gradient
% is above threshold
signal_filtered2 = signal_filtered;
signal_filtered2(TF) = 0;
figure(10),
plot(time,signal_filtered,'r',time,signal_filtered2,'b')
sound(signal_filtered2,Fs)

you will notice that the large spikes get's replaced by zeroes and so you will not be annoyed anymore with those klicks (red portion of the signal)

Jacob 2023-12-1

Interesting yeah I heard the extra noise in it but I was going to roll with it because I could actually hear the audio and not the tone in it! But thank you I will use this! Thank you for your help!