trouble making a notchfilter to remove a tone out of a wavefile
    조회 수: 5 (최근 30일)
  
       이전 댓글 표시
    
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일
				audiowrite('processed_SunshineSquare.wav', y_processed, Fs_processed);
soundsc(y, Fs);
None of those variables is defined in your code. Not even Fs.
채택된 답변
  Mathieu NOE
      
 2023년 11월 29일
        
      편집: Mathieu NOE
      
 2023년 11월 29일
  
      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
추가 답변 (0개)
참고 항목
카테고리
				Help Center 및 File Exchange에서 Get Started with Signal Processing Toolbox에 대해 자세히 알아보기
			
	Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!





