44.1 kHz Sample Rate

68 views (last 30 days)
Edy Nastase
Edy Nastase on 22 Jul 2022
Commented: Star Strider on 22 Jul 2022
Hello,
I have a very brief question, which bothers me a lot. So, I was playing around with dsp.SineWave and dsp.FIRFilter/dsp.LowPassFilter, but whenever I introduce the signal into the filter, the magnitude peak is still present, just at a lower dB.
After some short analysis, I realized that the "sampling rate" must be 44.1 kHz, and any other sampling rate will result in the peak returning.
Could anyone please explain this phenomena regarding sinusoidal waves and filtering?

Answers (2)

Star Strider
Star Strider on 22 Jul 2022
Could anyone please explain this phenomena regarding sinusoidal waves and filtering?
Digital filters are designed to work with the same sampling frequency as the signals they are filtering. Using them with signals sampled at any other sampling frequency will not produce the desired results. Apparently the filter you are using was designed with a sampling frequency of 44.1 kHz, so the signals you are using with it must have that same sampling frequency.
  5 Comments
Star Strider
Star Strider on 22 Jul 2022
%% Specify some basic parameters
N = 120; % filter order
Fs = 48e3; % sampling frequency
Fp = 8e3; % passing band
Ap = 0.01; % attenuation in the passband, in dB
Ast = 80; % attenuation in the stopband, in dB
Fst = 10e3; % Stopband frequency
% calculate the ripples in linear units (convert from dB)
Rp = (10^(Ap/20) - 1)/(10^(Ast/20) + 1); % passband ripple in linear units
Rst = 10^(-Ast/20); % stopband ripple in linear units
% 'minorder' = indicates that the lowest filter order
% f = frequency range ( [0 1], where 1 is Nyquist)
% a = amplitude of the signal [0 0 1 1] would give back a high pass
% DEV = deviation in the pass/stop band ripples
B = firgr('minorder', [0 Fp/(Fs/2) Fst/(Fs/2) 1], [1 1 0 0], [Rp Rst])
B = 1×157
1.0e+00 * -0.0000 0.0000 -0.0000 -0.0000 0.0000 0.0000 -0.0000 -0.0000 -0.0000 0.0000 0.0001 -0.0000 -0.0001 -0.0000 0.0001 0.0001 -0.0001 -0.0002 -0.0001 0.0002 0.0003 -0.0001 -0.0004 -0.0002 0.0004 0.0005 -0.0000 -0.0007 -0.0005 0.0005
%% Filter Data %%
f1 = 5000;
f2 = 15000;
sine = dsp.SineWave('Frequency', f1, 'SampleRate', Fs, ...
'SamplesPerFrame', 4000);
sine_1 = dsp.SineWave('Frequency', f2, 'SampleRate', Fs, ...
'SamplesPerFrame', 4000);
% Form the FIR Filter using the filter coefficients
% dsp.FIRFilter can be used instead of "filter"
% It has the advantage of managing the state
lowpass_filter_fir = dsp.FIRFilter('Numerator', B)
lowpass_filter_fir =
dsp.FIRFilter with properties: Structure: 'Direct form' NumeratorSource: 'Property' Numerator: [-7.6044e-06 2.7888e-05 -2.8541e-05 -1.0648e-05 2.0291e-05 2.4305e-05 -1.0894e-05 -4.1386e-05 -1.3330e-05 4.9489e-05 5.4346e-05 -3.0304e-05 -9.8008e-05 -3.0322e-05 1.1195e-04 1.2487e-04 -5.7249e-05 -2.1042e-04 -8.3246e-05 … ] InitialConditions: 0 Show all properties
figure
freqz(B,1,2^16,Fs) % 'firgr' Filter
figure
freqz(lowpass_filter_fir.Numerator,1,2^16,Fs) % 'lowpass_filter_fir' Filter
% return
% Displays the frequency-domain signals and
% the frequency spectrum of time-domain signals
% % spec_a = spectrumAnalyzer("SampleRate", Fs, "PlotAsTwoSidedSpectrum", false);
% % spec_a.ChannelNames = {"Input", "Output"};
% for i = 1:1000 % 1000 frames analyzed
% x = sine() + sine_1() + 0.5*randn(sine.SamplesPerFrame, 1); % add noise
% y = lowpass_filter_fir(x); % low-pass filter the signal
% spec_a(x, y); % use spectrum to showcase the signal
% end
% release(spec_a) % release the value in real time
Large parts of the DSP Toolbox code will not run online because the spectrumAnalyzer funcion will not work online.
However, both of the filters have cutoff frequencies of 10,000 Hz, and the sine waves have frequencies of 5000 and 15000 Hz. The 5000 Hz signal will pass, and the 15000 Hz signal will be filtered out —
t = linspace(0, 1E+5, 1E+5)/Fs;
sine = sin(2*pi*t*5000);
sine_1 = sin(2*pi*t*15000);
x = sine + sine_1;
x_filt = filtfilt(B,1,x);
figure
subplot(2,1,1)
plot(t, x)
xlim([0 0.005])
grid
title('Original Signal')
subplot(2,1,2)
plot(t, x_filt)
xlim([0 0.005])
grid
title('Filtered Signal')
Fn = Fs/2;
L = numel(t);
NFFT = 2^nextpow2(L);
FTx_xf = fft([x(:) x_filt(:)],NFFT)/L;
Fv = linspace(0, 1, NFFT/2+1)*Fn;
Iv = 1:numel(Fv);
figure
subplot(2,1,1)
plot(Fv, abs(FTx_xf(Iv,1))*2)
grid
title('Fourier Transform Of The Original Signal')
subplot(2,1,2)
plot(Fv, abs(FTx_xf(Iv,2))*2)
grid
title('Fourier Transform Of The Filtered Signal')
The specified filters work as desired!
.

Sign in to comment.


Paul
Paul on 22 Jul 2022
Edited: Paul on 22 Jul 2022
The main question is: why is there a magnitude peak, even after the signal has been filtered?
If the input to a stable, linear time invariant filter is a sine wave, the steady state ouput will be sine wave of the same frequency as the input. The steady state output amplitude and phase are determined by the gain and phase of the filter at the frquency of the input. So unless the filter gain is exactly zero at the input frequency, the steady state output will be at the same frequency. If if much attenuated, we'd still expect a peak in the output at the input frequency. Because the filter is linear, if the input is the sum of two sine waves, the output will be the sum of the outputs due to the individual inputs.
Looking at the code
clear all
close all
clc
%% Specify some basic parameters
N = 120; % filter order
Fs = 48e3; % sampling frequency
Fp = 8e3; % passing band
Ap = 0.01; % attenuation in the passband, in dB
Ast = 80; % attenuation in the stopband, in dB
Fst = 10e3; % Stopband frequency
% calculate the ripples in linear units (convert from dB)
Rp = (10^(Ap/20) - 1)/(10^(Ast/20) + 1); % passband ripple in linear units
Rst = 10^(-Ast/20); % stopband ripple in linear units
% 'minorder' = indicates that the lowest filter order
% f = frequency range ( [0 1], where 1 is Nyquist)
% a = amplitude of the signal [0 0 1 1] would give back a high pass
% DEV = deviation in the pass/stop band ripples
B = firgr('minorder', [0 Fp/(Fs/2) Fst/(Fs/2) 1], [1 1 0 0], [Rp Rst]);
%% Filter Data %%
f1 = 5000;
f2 = 15000;
sine = dsp.SineWave('Frequency', f1, 'SampleRate', Fs, ...
'SamplesPerFrame', 4000);
sine_1 = dsp.SineWave('Frequency', f2, 'SampleRate', Fs, ...
'SamplesPerFrame', 4000);
% Form the FIR Filter using the filter coefficients
% dsp.FIRFilter can be used instead of "filter"
% It has the advantage of managing the state
lowpass_filter_fir = dsp.FIRFilter('Numerator', B);
The frequency response of the filter at the two frequencies of interest is
h = freqz(lowpass_filter_fir,[f1 f2]/Fs*2*pi)
h =
0.7071 - 0.7071i -0.0001 - 0.0001i
The gain in dB
db(abs(h))
ans = 1×2
-0.0000 -80.7467
We we see that sine and sine_1 should be attenuated by 0 and 80 dB respectively, which is what is seen in spec_a (which can't be run on Answers apparently).
% Displays the frequency-domain signals and
% the frequency spectrum of time-domain signals
% spec_a = spectrumAnalyzer("SampleRate", Fs, "PlotAsTwoSidedSpectrum", false);
% spec_a.ChannelNames = {"Input", "Output"};
% for i = 1:1 % 1000 frames analyzed
% x = sine() + sine_1() + 0*0.5*randn(sine.SamplesPerFrame, 1); % add noise
% y = lowpass_filter_fir(x); % low-pass filter the signal
% spec_a(x, y); % use spectrum to showcase the signal
% end
% release(spec_a) % release the value in real time

Community Treasure Hunt

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

Start Hunting!