Choosing right cut-off frequency for a LP filter in upsampler I'm implementing a upsample function in Matlab but it's not perfect right now,for reasons I'm not sure. Here is my code:

U=5; %upsampling rate
N = U*length(x);

% Init
x_up = zeros(1,N);

% Zero-pad
ind=1;
for n=1:N
    if mod(n,U)==0
       x_up(n)=x(ind);
       ind = ind+1;
    end
end

%% Design a LP filter
tol = (1/U)/2;
[n,fo,mo,w] = firpmord([1/(2*U)-tol 1/(2*U)+tol],[1 0],[0.01,0.01]);
 b = firpm(n,fo,mo,w);

x_up = [zeros(1,n),x_up]; %To avoid group-delay
x_up = conv(b,x_up,'full');
x_up(abs(x_up)<0.01)=0;
% Remove leading and tailing zeros
x_up = x_up(find(x_up,1,'first'):find(x_up,1,'last'));
x_up = x_up(1:N);

And these are the results:

From theory perspective, everything looks like fine to me but how can I make my upsampler perfect?

View your frequency response after your low pass filter on a dB scale to better show the limitations of your filter.

Use a multiband filter with the least squares algorithm for an optimized rejection filter for zero-fill interpolation. This will concentrate the rejection to be specifically where the images are that need to be removed.

Given your original signal extends to the band edge, you will have significant challenges interpolating this with zero insert and filtering since the filter complexity is driven by the width of your transition band. I give details below on such an interpolator filter design - notice specifically what defines the signal you want to pass and the signal you want to reject. To simplify that and make this approach feasible, consider first low pass filtering your waveform to reduce the high frequency components while still being acceptable for what you want to have in the time domain. Then you can interpolate using the zero-insert and multi-band filter design approach. The low pass step is not necessary (as the interpolate filter can be designed to do the same thing) but helpful as you will be able to observe your signal prior to further processing to access if the time domain is still acceptable.

To the extent you can pass your signal within the pass band without any distortion, and completely eliminate the images, you will achieve perfect interpolation (this is not possible without an infinitely long filter and infinite delay so you can only approximate this-- given enough taps you can do very well).

This first slide shows how images are created through zero-insert. The digital spectrum is already periodic, zero insert just extends the sampling rate without changing the periodicity or original shape of the spectrum-- so when we interpolate by 4 for example as done here, the images that were originally around $F_s$, $2F_s$ become part of our spectrum in the first Nyquist zone at $F_s/4$ and $F_s/2$. These must be filtered out to complete the interpolation.

An optimized interpolation filter need only reject the images while not distort the passband signal of interest. A least squares (firls in MATLAB/Octave) FIR filter is an excellent choice and support multiband filters where we can concentrate the rejection where it is needed most, providing greater rejection for the same number of taps.


What language or framework would you use to develop and test natively algorithms meant to be ported on an embedded hardware later? I am trying to find a nice clean solution to this problem since some days now, and couldn't make myself to stick to a solution.

I need to develop some stream and block processing for audio, including cascaded biquad, FT and iFT, FIR, dynamic correction. The whole audio chain will be embedded either in a DSP or on a ARM running Xenomai Linux.

What I am looking for is a way to quickly compile and test bits of codes, getting audio in and out, and in a language that will allow me to target a broad range of devices. (Would be C or Cpp surely).

I don't want to choose the definitive platform now because I don't know yet how much processing power will be required to make my full project runs.

Here are some tought:

I am using Gnu Octave but afaik it will not compile a real-time pipeline and the code is hardly portable.
I tried Audio Weaver from DSP Concepts (which is great) but the targets DSP are limited and the engine is proprietary.
I looked at Faust which is a very promising framework for audio (even able to generate Cpp code), but unfortunately it was designed with stream processing in mind and not meant to process blocks of data
I looked at some LabView + FPGA solution but the price and licensing is just WOW !
I looked at some IDE from Ti and AD but there also the price is too high.
Using a simple VST or ALSA framework might be ok but I don't really know where to start.

I finally settled on python, and here is why:

I had no reason to limit myself to a specific language or framework. The target plateform was not specified, but I needed to have proof of concepts for the project, and quick.

I choose python because of its versatility, the profusion of libraries, the speed over other scientific calculation packages, and the scientific community.

And so far so good.


FIR Filter phase shift and stability I designed an FIR highpass filter which is causing a phase shift to the output. The output signal is also not stable (?) on start.
Are these normal or am I doing something wrong?
If these are normal, what causes these?
(Attached is a diagram showing the input in blue and the output in red).
Thanks.

This is perfectly normal. At the beginning you see the transient effect because the signal suddenly starts and the filter was at rest before (zero initial condition). A causal filter will always add some phase shift to your signal. This phase shift is usually frequency dependent, but for a linear phase FIR filter the resulting delay is independent of frequency.

If $H(e^{j\omega})$ is the frequency response of the filter and $\phi(\omega)=\arg\{H(e^{j\omega})\}$ is its phase response, then the phase shift experienced by a sinusoid with frequency $\omega_0$ is simply $\phi(\omega_0)$, i.e. the corresponding time delay equals $-\phi(\omega_0)/\omega_0$.


Fundamental frequency Can someone help me:

I have a signal which consist of three frequencies; for example, $f_1=800\,\text{Hz}$, $f_2=1050\,\text{Hz}$, $f_3=1600\,\text{Hz}$. How can i compute the fundamental frequency?

Also, is there a method called "frequency histogram method" which is used to compute fundamental frequency? Can someone explain it to me?

Thanks.

A signal which is the sum of sinusoids with different frequencies $f_i$ is only periodic if

$$f_i=k_if_0,\quad i=1,2,\ldots$$

is satisfied for some $f_0$ and integer values $k_i$. Or, in other words, the frequencies $f_i$ must be rational multiples of each other. In your example, the signal is periodic but its fundamental frequency is not $f_1=800\,\text{Hz}$ because $f_2\neq k_2f_1$ for integer $k_2$. The fundamental frequency is $f_0=50\,\text{Hz}$, even though this frequency does not explicitly occur in the signal. However, note that its period equals $1/f_0$.


How Is Mixed Norm ($ {L}_{1, 2 }$) Better than $ {L}_{1} $ Norm for Sparse Representation? Using $ {l}_{1} $-norm regularization for the purpose of achieving sparsity of the solution has been successfully applied a lot. But many times, I found the paper using mixed-norm instead of $l_1$-norm.

Considering the mixed norm $ {L}_{p,q}$ norm defined as:

$$ {\left\| a \right\|}_{p,q} = \left( \sum_{i} \left( \sum_{j} {\left| {a}_{i,j} \right|}^{p} \right)^{
\frac{q}{p} } \right)^{ \frac{1}{q} } $$

How is mixed-norm better than $l1$-norm for sparse representation? Most of what I have seen is $l_{1,2}$.

Reference, for example,

Sparse Nonnegative Matrix Factorization for Clustering eq.12
Compressive Sensing of Jointly Sparse Signals as a Method for
Dimensionality reduction of Mass Spectrometry Data eq.7
L1/2 Sparsity Constrained Nonnegative Matrix Factorization for Hyperspectral Unmixing

The mixed norm allows you to impose some simple structure in the solution matrix. Using your example with $p=q=1$ then this means the solution could have arbitrary elements set to non-zero coefficients. This would not impose any structure in the solution. In this case the $L_{11} $ norm just sums the absolute value of all the matrix elements - similar to the Frobenius norm.

As an example, consider that $L_{11}$ norm of the following matrices are the same:

$\begin{bmatrix} 0&2 \\ 0 & 0\end{bmatrix}, \quad \begin{bmatrix} 0&0 \\ 2 & 0\end{bmatrix}, \quad \begin{bmatrix} 2&0 \\ 0 & 0\end{bmatrix}, \quad \begin{bmatrix} 0&0 \\ 0 & 2\end{bmatrix}$.

So there is no preference as to which element is set to non-zero

Instead, if $p=2$ and $q=1$, then the objective function is minimizing the $l_2$ norm of the columns and then the $l_1$ norm over the that. This means you're looking for an $a$ matrix with a sparse number of column vectors and where each column vector has a small $l_2$ norm. Of course this assumes you are using the $l_1$ norm as a proxy for sparsity.

Consider the $L_{2,1}$ and $L_{1,2}$ norms of following matrices:

$\begin{bmatrix} 2&0 \\ 2 & 0\end{bmatrix}, \quad L_{2,1} = \sqrt{8} \approx 2.8, \quad L_{1,2} = 4
$

and

$\begin{bmatrix} 2 & 0\\ 0 & 2 \end{bmatrix}, \quad L_{2,1}= 4, \quad L_{1,2} = \sqrt{8}
$

So minimizing the $L_{2,1}$ norm tends to prefer the structure of the first matrix over the second one, while minimizing the $L_{1,2}$ norm tends to prefer the structure second matrix over the first.

Thus the mixed norm allows you impose a soft constraint on the structure of the sparse solution that you're looking for. It's relatively easy to see that be you could also set it up to look for a sparse set of row vectors, each with a small $l_2$ norm. Deciding which to use is really dependent on how you set up your equations.


Algorithm that Fills the White Holes Within the Character Image I have this degraded image :

Using paintbucket tool in Paint.NET Windows application, I fill the character image black.

My questions are :

What algorithm can we use to fill up the holes within the character image?
Is there any corresponding built-in function in Python that fills up the white holes in a character image?
What probable algorithm does the paint bucket tool in Paint.Net uses ?

Thanks.

Based on the blog post - The Paint Bucket in Paint.Net 4.0 (Video) I can tell it uses some edge detection to handle similar colors within a piece wise smooth area.
More information is given in the Paint Bucket Tool documentation.
Usually the way it can be implemented is by defining color metric. How far a color is form another color. If it within the tolerance, it is within the mask.

Morphological operation can fill small holes.
I am not sure, but it won't be too hard to implement something similar.
Probably they convert the RGB data to HSL / HSV Color Model. Then just create a mask based on the distance between the reference point to the surrounding colors.  Now the mask creates a logical (Black and White) mask to be filled.
In case of Paint.NET the algorithm should be something like Flood Fill.

So the steps are:

Convert image to the Color Model of choice (HSV / HSL or even LAB).
Generate a mask using color similarity (Distance Function) on the Color Model. The distance function is parameterized by some Threshold parameter.
Use Flood Fill like algorithm to fill the pixels in the mask using the target color. The mask is basically the stopping condition for the Flood Fill operation.

In MATLAB you have imfill() and Flood Fill.
For something more advanced you may use adaptive morphological functions like Active Contours and Morphological Snakes.


Is DCT basicly a filter? I am reading about audio processing for my university project and I found out that DCT basicly filters out some frequencies. I wanted to use for example high-pass filter and then DCT for compression, but isn't dct a filter?

What's the difference?

DCT and DFT can be thought of as a bank of bandpass filters.  So it's sort of a filter, but it doesn't remove any part of the signal.  The whole signal is still there, just in a different format.

Yes, you could do the DCT on the signal and then high-pass filter by multiplying by the DCT of a high-pass filter's impulse response.  I'd suggest using DFT/FFT first though; it's a little more intuitive, though for a final application DCT has better energy-compaction, and so can compress signals better.


Small shift in frequency in FFT I'm trying to compare the estimation of the power spectral density (in this case of a fake cosine signal) by two different methods: 1st, MATLAB's fft function, and 2nd, a direct calculation from the definition, $\int f(t)\cdot e^{-i\omega t}\mathrm{d}t$. The problem I find is that, for whatever reason, the FFT result is slightly shifted towards the right in frequency, and I'm not able to see why.

I'm adding the full code so that if anyone wants to try it, it should directly work and plot the results in MATLAB, something like

>> compare_FT(10);

As is clear from the code, the peak should be exactly at $\omega = 6$: that is indeed the case for the calculation from definition, but not exactly so for the FFT.

function[] = compare_FT(fs)

%% signal parameters (now fake, then from lab)
f = 6/(2*pi);
w_d = 2*pi*f;
T = 100;
num_samples = T*fs;

%% vector definitions
t = linspace(0, T, num_samples);
xx = cos(w_d*t); % the "fake" signal
dt = t(2) - t(1); % or, equivalently, dt = 1/fs
% omega vector for DTFT in a specific range:
ww_size = num_samples; % the bigger this size, the more points in the estimated FT
w_start = 5; % start calculating the FT at this omega
w_end = 7;  % stop at this omega
ww = linspace(w_start, w_end, ww_size);
Sx = zeros(1, ww_size); % prealocate space for integration

% DTFT method
for i = 1:ww_size
    int_value = trapz(xx.*exp(-1i*ww(i)*t))*dt;
    Sx(i) = int_value.*conj(int_value)/T; % estimate PSD
end

%% FFT method
n_fft = 2^14; % number of points in FFT (more points -> better FT discretization)
signal_fft = fft(xx, n_fft); % the fft itself
X2_fft = abs(signal_fft).^2; %signal_fft.*conj(signal_fft);
L = length(X2_fft);
Sx_fft = X2_fft(1:L/2+1)/(fs*num_samples); % proper scaling FFT to estimate PSD
w_fft = 2*pi*fs*(0:L/2)/L;

%% Plot results to compare both methods
figure(3);
clf;
box;
hold on;
% Plot theoretical value at top
plot(ww, ones(1,length(ww))*T/4, '--', 'Color', 'g', 'LineWidth', 3.5);
% Plot "direct" calculation
plot(ww, Sx, '.', 'Color', 'b', 'LineWidth', 3.5);
% Plot "FFT" calculation
plot(w_fft, Sx_fft, '-', 'Color', 'r', 'LineWidth', 3.5);
set(gca,'FontSize',18,'FontName', 'CMU Sans Serif');
xlabel('\omega', 'fontsize', 24, 'FontName', 'CMU Sans Serif');
ylabel('$S_x$', 'fontsize', 24, 'Interpreter', 'latex', 'FontName', 'CMU Sans Serif');
title('DTFT');
xlim([w_start, w_end]);
grid on;
end

Your vector

t = linspace(0, T, num_samples);

is the problem. It should be

t = (0:num_samples-1)/fs;

Note that, even though you claim it in a comment in the code, in the current version, dt = t(2) - t(1); is not equal to 1/f_s.


Applying the CUSUM algorithm to a correlated random process As far as I know, the CUSUM algorithm is meant for detecting change points on discrete-time uncorrelated random processes.

For instance, to apply the CUSUM algorithm to a discrete Gaussian process, we must know for sure that each sample is statistically independent from the others. I have seen this assumption in the CUSUM algorithm demonstration.

However, I have not found the CUSUM alternative for a discrete-time correlated random process. Let's say, for a realization of a band-pass filtered Gaussian process $X[n]$ with a changing variance at a certain sample $n_c$.

What to do in such a scenario? Can I still apply CUSUM?

If you look at the derivation of CUSUM from first principles, the idea is to construct a likelihood ratio test statistic and wait till it exceeds a threshold. In case of CUSUM, the assumption is that samples are uncorrelated and so the test statistic simplifies to a sum. Since you are dealing with a different statistical change model, you can try writing a likelihood ratio statistic and see if you can simplify it to anything useful. There are also other algorithms out there for testing "shifts in variance" or "shifts in spectral content" for ARMA processes which you may want to read about.

See also: Ch 6 and 7 Detection of abrupt changes: theory and application
M Basseville, IV Nikiforov - 1993.


In 1D DCT, why is the input a vector? The question is specific to this document: Image Compression.

It is chapter 6 from book "A First Course in Applied Mathematics" by Jorge Rebaza. It is, a to the point explaination of DCT that is easier for an engineer to grasp. However, some things are not explained and thus here I am.

The questions are about the notation used when the book starts talking about 1D DCT; on pg.3 of the pdf it shows $x$ that is the input to be transformed. This $x$ is a vector i.e a matrix with $n$ rows and $1$ column. We then multiply this with $C$ that is $n \times n$ matrix to get a $n \times 1$ matrix $y$, that contains the DCT coefficients.

Q: Why is $x$ described fundamentally as a column matrix and not a row matrix on pg.3 of pdf? If the $x$ was row matrix i.e it was not transposed, making it $1 x n$ matrix, how would the equation to calculate $y$ change? Right now it is $y=Cx$.

Q: Why is the equation y=Cx and not written as $y'=Cx'$ where the $'$ symbol means transpose? Since we are describing both $x$ and $y$ as being transposes?

Q: When we move to 2D DCT (pg.8 in pdf) the equation is derived like this $y=C \cdot (Cx')'$ which gives $y=CxC'$. How come here we explicitly write the transpose of the operand $x$ but not in the 1D equation?

The descriptions are looking inconsistent.

I do find this chapter quite well-written. The row or the column visions of a vector are mostly matters of convention and convenience. One or the other is perfectly fine.   In signal processing, when using algebra, the column version (a vertical vector) is used quite often. If you multiply a col-vector $x$ by matrix $A$, you get a col-vector $y$:

$$ y = Ax$$

You can get the row-vector version by transposing:

$$ y' = x'A'$$

If you decide to work with row-vectors, you could as well work with $v=uB$, but the matrix on the left seems more easy to read (at least with western languages).

When you deal with an image or image block $X$, you can apply a transform on both the columns, at $Y_1=AX$. $Y_1$ is a matrix whose columns are the transformation of each column of $X$. Then, you can work on the rows of $Y_1$, by applying the formula  $Y_1 A'$: every row is now affect.ed by the transform. The final result is thus $(AX)A'$, or $A((AX')') = A(XA')$


How to Find Performance Ratio (PR)? If we have one image from BSD (The Berkeley Segmentation Dataset and Benchmark) and for this images Ground truth, how to calculate PR?

PR is the ratio of true to false edges and can be calculated by the formula:

  PR = True Edges (Edge pixels identified as Edges)/False Edges (Non
  edge pixels identified as edges) + (Edge pixels identified as Non-Edge
  pixels)

How to find those arguments in formula?

When you segment an image you have boundaries between segments.
Those boundary pixels (Just between 2 different segments) are the edge pixels above.


Resources for learning about filters from scratch? Topic of filters is very important in signal processing
So i want to learn about filters from scratch and then eventually reaching FIR and IIR filters
I have some idea about low pass,high pass,band pass,band stop types but often get confused how they relate to butterworth, chebyshev and other such types?
Apart from wikipedia, please recommend books chapters or websites for learning about filters

The best book I know that exclusively treats digital filters is Digital Filter Design by Parks and Burrus.
But any good textbook on signal processing has one or more chapters about properties and design of digital filters. And since the most common IIR filter designs are based on analog filters, these books usually also treat the classic analog filter approximations like Butterworth, Chebyshev, Cauer (elliptic), and Bessel.
Two good undergraduate texts are Digital Signal Processing by Proakis and Manolakis and Digital Signal Processing Using Matlab by Ingle and Proakis.
Also take a look at the answers to this very related question.


All-pass phase compensation only in a specific frequency range I am trying to design an all-pass filter in order to compensate the phase distortion at the output of an equalizer. To do so, I take the phase vector of the output, compute the real and imaginary part using  $e^{j\phi}$, then I do an IFFT and save the time-reversed output. The saved impulse response of the all-pass filter looks like this:

Problems start when I try to restrict the working range of the all-pass equalizer: for example, let's say that the phase compensation needs to work only in 400-16k Hz range, leaving the remaining frequencies uncompensated. To do so, I'm essentially zeroing the phase outside the working range before computing the real and imaginary parts. Probably this is not the correct way, since the resulting impulse response has several ringings at the start and the end of the plot, as you can see from the following picture.

Could you please help me to understand which is the problem here?
EDIT in response to @Dan Boschen's answer
This is the solution I implemented in Octave following your suggestions:
h = read_dat("heqm.dat"); % the IR to equalize
NFFT = length(h);

H = fft(h,NFFT);
ph = unwrap(angle(H));

x = 1:NFFT;
p = polyfit(x,ph,1);
f = polyval(p,x);
f = f-f(1);
figure; plot(ph); hold on; plot(f);

The corresponding plot is:

EDIT 2
The real(ifft(exp(1j*ph))) is:

Next, the definition of the mask:
f1 = 100;
f2 = 300;
half_mask = zeros(1,NFFT/2+1);
half_mask(f1:f2) = 1;
mask = [half_mask, fliplr(half_mask(2:end-1))];
idx = find(mask);

res = f;
res(idx) = ph(idx);
figure; plot(res);

which corresponds to:

Finally, the IR computation:
imp = real(ifft(exp(1j*res)));
figure; plot(imp)

which unfortunately has aliasing:

Did I do something wrong?

The following demonstrates an approach close to what the OP describes that does not lead to the time domain aliasing experienced. Specifically, the average delay of the filter is determined and for that linear slope of phase versus frequency- all frequencies that are not to be modified by the equalizer are to be replaced with that linear phase slope.  (Update: after reviewing OP’s specific waveform the time domain aliasing still exists following this process, and I am concluding that the time span of the equalizer is not supporting the overall time delay given by the waveform  - see bottom of post for that detail)
Consider a phase response in the frequency domain as follows from an example 5000 bin DFT, together with a best fin linear line consistent with the average delay of the signal:

Below is the IFFT prior to any selection
imp = fft.ifft(np.exp(1j*phase))
plt.plot(fft.fftshift(imp))

We note that bin 0 represents "DC", bins 1 through 2500 the positive frequencies, and bins 2501 through 4999 represent the negative frequencies. (With bin 4999 the negative of bin 1, bin 4998 the negative of bin 2, etc).
If we consider this phase response to be the compensation we want, but only isolated to selective frequencies, we have to consider if the signal that we are compensating is real (after compensation, in which case the phase will be conjugate symmetric after compensation) in which case we need to include both positive and negative frequencies, or if the signal is complex and if we only wish to compensate the positive or negative frequency component.
Showing both cases to show that the approach handles both cases without aliasing issues:
First for a real result after equalization:
f1 = 800
f2 = 1500
N = 5000
mask = np.zeros(N)
mask[f1:f2+1] = 1
mask[N-f1:N+1 -f2] = 1
result = np.where(mask, phase, plin )

With the following result for the impulse response:

And again for the case of selecting only the positive frequencies:
mask = np.zeros(N)
mask[f1:f2] = 1
result2 = np.where(mask, phase, plin )

I do not know the specifics of what the OP had done, but the time domain aliasing shown is consistent with zeroing every other bin in the frequency domain. For example for the phase plot shown in the beginning if we zeroed every other sample a similar time domain aliasing would occur:
phase2 = np.copy(phase)
phase2[::2]= 0
imp4 = fft.ifft(np.exp(1j*phase2))
plt.figure()
plt.plot(fft.fftshift(imp4))

UPDATE:
I reviewed the OP's data and conclude from that the overall delay is exceeding the time duration of the samples given causing aliasing in the time domain to occur.
Observe the following results:
Out of personal preference, I include fftshift in my plots of the impulse response to center the plot at t=0, plotting positive and negative time as an equivalence to the processing in the frequency domain with positive and negative frequencies. Similar to getting aliasing for frequencies above $f_s/2$ (negative frequencies) we can have similar effects in the time domain for delays above $T/2$ where $T$ is the total time duration of the signal being processed.
To illustrate, below is the process to plot the impulse response for a vector starting at t=0 with sample 1 followed by 1023 zeros (matching length of the OP's example, and coded on Octave (and instead of taking the real, I plot the complex magnitude):
xaxis = -512:511
plot(xaxis, fftshift(abs(ifft(exp(j*zeros(1024,1))))))):

Same processing for impulse response of OP's signal and it is clear that the overall delay has been further shifted beyond half the total time duration which would result in aliasing (consider the impulse response of the equalizer would ideally be centered in the span of the filter while in this case it has gone beyond the filter edge):
plot(xaxis, fftshift(abs(ifft(exp(j*phi)))))):

In general, for every delay of one sample, the phase will vary linearly from 0 to $-2\pi$ from DC to the sampling rate. The OP's phase is varying approximately -450 radians over 1024 frequency bins. $450/(2\pi)$ corresponds to a delay of 72 samples
plot(xaxis, fftshift(abs(ifft(exp(j*zeros(1024,1))))))):

Here is the OP's original waveform along with a 72 sample delay (which would then wrap cleanly at the positive/negative frequency boundary as well, but still resulted in the time domain aliasing the OP has computed):

Here is the impulse response computed from the phase of the 72 sample delay to confirm the prior processing:
plot(xaxis, fftshift(abs(ifft(exp(j*f)))))

Conclusion: I am rather surprised that the waveform would have that much additional delay due to the variation of the phase from the best linear fit (as an analogy to frequency translation occurring (in the frequency domain) from phase rotation in time, here we clearly see evidence of time translation (in the time domain) from the same phase rotation process apparently occurring in the frequency domain. Given this, the impulse response is being translated further than the original time sample can support, and the solution would therefore require a longer time sample which should then eliminate the aliasing. (Sometimes it is easier to reverse time and frequency to understand what may be occurring; as it may be more intuitive for many to understand how this could occur with frequency domain aliasing from undersampling--- here it is the same thing except instead of increasing the sampling rate we need to increase the time duration).


How to design a filter that can filter out noise accurately, after setting the parameters of the filter using standard signal? I meet a problem with designing a filter. I have two different instruments that could measure the same AC signal (usually ~200hz, always <1kHz), A and B.  A can carry out signal measurement during the normal operation of the instrument.
B can only be used during instrument calibration, and A can also be used at this time. B can accurately measure the signal. The measurement accuracy of A is lower than that of B.

I want to design a real-time filter for A to filter out the interference noise. The parameters of the filter are adjustable. When the instrument is calibrated, both A and B are used to measure the signal.
The measurement result of B is taken as the standard result to calculate the parameters of the filter adaptively. After setting the parameters of the filter, A is used to measure the signal without B.

I have investigated a variety of filters and found that the adaptive filter is quite suitable for my needs. However, the adaptive filter needs to input two signals, while I can only provide the measurement result of A during the operation of the instrument. So the adaptive filter is not completely suitable for my needs.

I don't know how to solve this problem. Is there any method that suits my needs?
Thank you very much!

If you can use a higher precision instrument to measure the signal, then, if the signal maintains it's characteristics, you can build a system's model and use a Kalman filter for the "not so good" input channel.


Find out when an audio file is actually recorded? Is there any technology to help find out when an audio file is actually recorded?

Here is the problem:

Given an audio file, how can we find out when it is recorded? File creation time is useless because it can be altered.

Baring information content that contains verifiable historical landmarks, you can't date data .  Anything can be altered in a file, including the recorded content, so the data itself is no more reliable than a time stamp in a file header.  The normal assumption about data files is that they can be faithfully reproduced over and over without degradation, so there is no time information in the data itself.  You may be able to infer rough time estimates about the age of a recording based on certain qualities in the recording (assuming that the data has been unaltered).  For example a digital reproduction of an analog recoding that was originally mastered on a particular media has qualities that may be discernible upon analysis.  If you can identify qualities that indicate the type of equipment used to capture the audio and type of media the audio was captured on, you may be able to infer that the recording is from a certain time period.


Example of Fourier Transform not existing for real-life signals? I got curious based on this question here, but basically, is there ever a real-life signal that exists where its Fourier transform does not exist? If a signal is not finite energy, then its Fourier Transform does not exist, so what might be an example, (if any), of such a signal in real-life?

All real-life signals are finite energy. The universe contains a fixed (and finite) quantity of energy, which has been unchanged since it came into being.

A signal's energy is given by

$E =\int_{-\infty}^{\infty}|x(t)|^2dt$

Thus, the only way to make a signal's energy go to infinity is to allow it to continue for infinite time or reach an infinite peak level. While useful in mathematical and/or physics theory, neither of these is possible in reality.


How to segment with false edge in in image? I try to resolve false edge problem in image segmentation. Let see its defination in Figure. The false edge problem often occurs in inhomogeneity image, especially for medical image. I am finding the method or feature to resolve it and achieve real edge. I try to apply some existing method example: localizing active contour, graph cut but they are not effective. I also apply gradient or entropy in here. Could you suggest some idea or feature or method to resolve it?Thank all . You can download test image at https://www.dropbox.com/sh/ttymp4vnpjhbvfx/2iacME5zYp . Note: orginal image is "input.mat". Its type is 16 bit

Let see example,In my marked image, B is real edge. We only analyse A point. At the A point, the intensity bigger than object region, but it not so big. Assume intensity of A point is 300, in which object intensity is about 260-280 (inhomogeneity property). So we can said A point is background, if we calssify it into two region. Sorry, you are true, the question must be update that A is true edge and C is false edge(because C is create by other region and background); true edge is created by object region and background

You could try "Non-maxima suppression" This is usually used on edge detectors to thin the edge to only one pixel. When there are lots of edges close together you basically pick the one with the largest gradient, and assume that is the true edge

Another technique that might be applicable is double thresholding. Some edge detectors (canny for instance) actually perform multiple detection's with different parameters.  "Strong edges" are found in the the higher threshold (and probably the lower one as well) weaker edges may only be found using the lower threshold, and are only "true" edges if they connect to one of the stronger edges. Maybe you could use your current edge detector as the weaker classifier, and then make one that is harder to satisfy. This might eliminate the false positives that are only detected by the weak detector (i hope that makes sense)

good luck


Which Transformation, Or Similarity Metric, Is Rotation, Shift and Scale Invariant? I have an algorithm to to detect copy-move forgeries in images. I used block matching to detect regions of an image that were forged with copy-move forgery and highlighted the alleged areas. For block matching, I took the DCT of each block and normalized with an extended 16x16 quantization matrix.

The problem:

this approach does not work if the copied region is rotated or scaled.

My question:

What transformation/similarity metric would be invariant to rotations, shifts, and scaling? (Alternatively, is there another method available for solving this problem?)
Additional points if the method is already available in MATLAB.

Thanks for the help.

The code for this project is available here: copy-move forgery detection.

Edit: I am aware of the Radial Gradient Transform, but I don't think it is invariant to all three.

Use a methods called Error Level Analysis (ELA). It can solve your problem as many forensic applications use it to detected doctored reigns in digital images.


how can I analyze a basic sine wave for its frequency phase and amplitude at a particular point in time in matlab How can I analyze a basic sine wave for its frequency phase and amplitude at a particular point in time in MATLAB? Are there any tools to do that? I am trying to do additive synthesis of inharmonic sounds in MATLAB, but having a great difficulty to extract the needed phase amplitude and frequency data of the partials I am trying to generate. Can someone shed some light regarding signal analysis in MATLAB?

If your the unknown signal $x(n)$ is modeled as:
$$x(n)=A \sin(2 \pi f n+ \phi)+w(n)$$
and you want to estimate $A$,$f$,and $\phi$ accurately, you can use least square estimation. Unfortunately the cost function is nonlinear. You can use nonlinear least square in MATLAB to find the parameters as follows:

Make a cost function:

f=@(A,f,phi) x[n]-A sin(2*pi*  f *n+ phi)

and use

p0=[A0,f0,phi0];

p = lsqnonlin(f,p0);

to find the unknown parameters.
Note that the optimizer will have a hard time finding $f$ as the problem is not convex. So it is best if you can give an initial estimate of frequency by using a method like fft.

If the frequency is known then the problem can be converted to a linear estimation as:
$$x(n)=A \sin(2 \pi f n+ \phi)+w(n)$$
$$=A \sin(2 \pi f n) \cos (\phi) +A \cos(2 \pi f n) \sin (\phi)+w(n)$$
$$=p_1 S[n] +p_2 C[n]+w(n),$$
where $p_1=A \cos(\phi)$ and $p_2=A \sin(\phi)$ are unknown parameters, and $S[n]$ and $C[n]$ are known.


Determine group delay from the frequency response of FIR filter I have a issue to find out, is there any way that one can determine a group delay from frequency response of the FIR filter. Let's say I have I type FIR filter with order of 51. Frequency response obtained from freqz MATLAB's function is presented below.

I know group delay is negative derivative of linear (in this example) phase response, but how can one determine that from the picture above.

The group delay is the negative derivative of the phase response as the OP has stated, and specifically for the delay of one clock sample, the phase will go linearly negative to $2\pi$ radians as the frequency goes from 0 to $f_s$
From the picture we see the phase is going approximately 800 degrees at a frequency of  $.17\pi$ rad/sample  (where $2\pi$ rad/sample is the sampling rate). So this would be equivalent to
$800/360*(2\pi)$ rad /$(.17\pi)$ rad/sample = $26$ sample delay
For linear phase filters, the group delay is half the order which would be 25.5 in this case (we can't resolve that from the graphic alone). The OP wrote that the order of the filter was 51 which implies 52 taps.


Finding the correct $x[n]$ from a given DFT $x[k]$ I'm struggling to solve this question from an exam. The question gives you a set of IDFTs from a given DFT, only one of those is the correct.

The first thing I did is notice that $2\sin(\pi k)$ will always be zero for every $k$ so I ignore it.

Then I proceed to do the IDFT with the rest as $x[n] = x[k]e^{2j(2 \pi /N)n}$

And there is where i get lost and confused cause i know that the cos can be expressed in the Euler form and also:

l -> lambda;   cos(wn) = l[n-m] + l[n+m] <- Modulation.

So, when i see the given answers... i think i am missing something.

Thanks for your time.

Let's solve the exercise. It's asking for your awareness of periodicity in DFT and some trigonometry.

First, beacuse of the question content, let's state the following :
$$\delta[n-d] \overset{DTFT} \longleftrightarrow e^{-j \omega d } \overset{DFT} \longleftrightarrow e^{-j \frac{2 \pi}{N} d k } $$

Then given our DFT $X[k]$ work it backwards into $x[n]$ as a sum of delayed impulses.

$$
\begin{align}
X[k] &= 2 \cos(\frac{\pi k}{2}) + 2 \sin(\frac{14 \pi k}{8}) + \sin(\pi k) \\
X[k] &= 2 \cos(\frac{\pi k}{2}) + 2 \sin( \left( 2\pi - \frac{2\pi}{8} \right) k) + 0 \\
X[k] &= 2 \cos(\frac{\pi k}{2}) + 2 \sin(- \frac{2\pi}{8} k) \\
X[k] &= 2 \cos(\frac{\pi k}{2})  - 2\sin(\frac{\pi}{4} k) \\
\end{align}
$$

We have used basic trigonometry to reach the above line. Now lets use Euler's formula to expand the cosines and sines into complex exponentials so that we can use the corresponence of the first line to deduce the time domain impulses.

$$
\begin{align}
X[k] &= 2 \cos(\frac{\pi k}{2})  - 2\sin(\frac{\pi}{4} k) \\
X[k] &= 2 \left( \frac{ e^{j \frac{\pi}{2} k} + e^{-j \frac{\pi}{2} k} }{2} \right)  -2 \left( \frac{ e^{j \frac{\pi}{4} k} - e^{-j \frac{\pi}{4} k} }{2j} \right) \\
X[k] &= \left( e^{j \frac{\pi}{2} k} + e^{-j \frac{\pi}{2} k}  \right)  +j \left( e^{j \frac{\pi}{4} k} - e^{-j \frac{\pi}{4} k} \right) \\
\end{align}
$$

Now the correspondence between $\delta[n-d]$ and its $N=8$ point DFT is:
$$\delta[n-d] \overset{DTFT} \longleftrightarrow e^{-j \omega d } \overset{DFT} \longleftrightarrow e^{-j \frac{2 \pi}{8} d k } = e^{-j \frac{\pi}{4} d k } $$

And by using ths relationship we can deduce from the last line that the impulses are:

$$
\begin{align}
X[k] &= e^{j \frac{\pi}{2} k} + e^{-j \frac{\pi}{2} k}  +j e^{j \frac{\pi}{4} k} - e^{-j \frac{\pi}{4} k} \\
x[n] &= \delta[n+2] + \delta[n-2] + j\delta[n+1] -j \delta[n-1]
\end{align}
$$

One last step remains to complete the answer; representation of $x[n]$ in the range $0 \leq n \leq 7$ as it is the valid range for the formal relationship between $x[n]$ and its $N$-point DFT $X[k]$. Then we should map those impulses outisde of this range into it. As you can see $\delta[n+1]$ and $\delta[n+2]$ are outside of this range and should be mapped into $0 \leq n \leq 7$ by shifting them right by $N=8$. This shift is the consequence of the inherent periodicity in both $x[n]$ and $X[k]$ as $x[n-N] = x[n]$ etc. So we decude that $x[n]$ in the range $0 \leq n \leq 7$ is:

$$x[n] = \delta[n+2-8] + \delta[n-2] + j\delta[n+1-8] -j \delta[n-1]$$
$$x[n] = \delta[n-6] + \delta[n-2] + j\delta[n-7] -j \delta[n-1]$$

which is shown in the option b.


maximum output of a LTI system within linear limits In book , 'fundamental of signal and system'by M.J ROBERT it is written that

  Since no practical system can ever produce an unbounded response, strictly speaking, all practical systems are stable. The ordinary operational meaning of BIBO instability is a system described
  approximately by linear equations that would develop an unbounded response to a bounded excitation if the
  system remained linear. Any practical system will become nonlinear when its response reaches some large
  magnitude and can never produce a truly unbounded response. So a nuclear weapon is a BIBO-unstable system
  in the ordinary sense but a BIBO-stable system in the strict sense. Its energy release is not unbounded even
  though it is extremely large compared to most other artificial systems on earth

For a given LTI system, what parameter of system along with input which  determines the maximum magnitude of response after which it becomes non linear and how to calculate this maximum magnitude of response mathematically?

For a given LTI system, …, after which it becomes non linear

never. If you're modeling a system to be LTI, it's LTI: Linear (and time-invariant).

You can of course model a system to be non-linear in general, but linear for a range  of input. But then the answer to your question is "that's in your own model, and we can't tell you"!


Signal processing using numpy python To process a .wav audio file with numpy (using fast Fourier transform algorithm). I want to process an audio signal at a particular interval with a sampling frequency 44100hz and sampling rate of 20ms and also want to store the data into a array by using numpy with fft of 1024

Can anyone help me with the code?

First of all, sampling frequency and sampling rate are synonymous. You mean sampling frequency of 44.1 kHz with a data length of 20 ms.

Second of all, 20 ms of data at 44.1 kHz will give you 882 points. Not enough for a 1024-point FFT.

You will either need to upsample to 51.2 kHz, that way 20 ms will give you 1024 points. Or you could append 142 zeroes to your data. Or, you could simply compute your FFT with the 882 points that you have. All 3 methods, will give you the same information. I personnaly recommend that you simply compute your FFT on 882 points.

As to the actual syntax, just google it, there are some pretty good examples.

edit : This seems like a good place to start

https://www.oreilly.com/library/view/elegant-scipy/9781491922927/ch04.html


Difference Between Convolution and Multiplication I read that multiplication is convolution in frequency domain. I also understand that convolution is just polynomial multiplication. Can somebody explain what are the advantages of doing convolution in frequency domain? Also, i would like to know how to multiply two signals in time domain and what is the significance of it.

Multiplication in your first sentence is term-by-term multiplication: $z[n] = x[n]y[n]$ for all $n$. Convolution, for discrete-time sequences, is equivalent
to polynomial multiplication which is not the same as the term-by-term
multiplication. Convolution also requires a lot more calculation: typically $N^2$
multiplications for sequences of length $N$ instead of the $N$ multiplications
of the term-by-term multiplication.

The key point of Fourier analysis is that term-by-term multiplication in one
domain is the same as convolution in the other domain.  So, in order
to calculate the results of a convolution, you can either do it directly,
using $N^2$ multiplications, or transform to the other domain, do a term-by-term
multiplication, and transform back. This requires two transformations
to go from one domain to the other, $N$ multiplications in the other domain,
and one inverse transformation to come back to the domain
where the convolution result is needed. This more complicated process
can, in fact, require less computation because the transformations can
be done very efficiently via the Fast Fourier Transform (FFT) algorithm
which requires about $N \log N$ multiplications. So, compare the computational
effort in calculating the three transforms and doing the term-by-term
multiplication
to $N^2$ to see if the FFT gives you a more efficient method of computing
a convolution.


DTFT of window function applied to input signal $$x[n] = cos(\omega_1n) + cos(\omega_2n)$$

$w[n] = 1/N$ for $0 \leq  n < N, 0$ for everything else

Find the DTFT of $y[n]=x[n]w[n]$ expressed by the DTFT of $w[n]$, $W(\omega)$

I was thinking that multiplication in time domain is convolution in frequency domain. But $X(\omega)$ is 0 except for two frequencies so the convolution integral is 0? Which can't be correct? What to do here?

It seems you are mixing discrete and continuous notations and it's unclear which it is for time and frequency. You use the term DTFT which would indicate both are discrete but you use the frequency index $\omega$ which would be continuous in the frequency domain.

Your time domain equation doesn't make any sense, since your left side depends on $n$ which doesn't show up in the right hand side at all. The way it's written $x[n]$ is simply a constant. I'm not just nit-picky, the exact way you write & define this can make a big difference.

  I was thinking that multiplication in time domain is convolution in
  frequency domain.

Correct

  But X(ω) is 0 except for two frequencies

Maybe. $X(\omega)$ is, but $X[k]$ may or may not be depending on the exact relationships of $N$, $\omega_1$, $\omega_2$ and the sample rate.

  so the convolution integral is 0?

Nope, it's not.

If the frequency domain is continuous  $X(\omega)$ is non-zero in a special way, so that the convolution integral is perfectly well defined and non-zero. If it is discrete, the convolution is a sum, not an integral.


wavelet packet transform and lifting scheme? so the lifting scheme is basically an alternative to performing the discrete wavelet transform with several advantages.

But here are three questions which I did not find an answer to:

is it possible to use lifting also for calculation of the wavelet packet transform or does it just apply to the normal wavelet transform?
is the wavelet packet transform also critically sampled like the dwt or is it redundant?
if using a level basis in the packet transform, i.e. a representation which uses only elements from the same decomposition level, the time-frequency plot is regularly discretized. How does this representation then differ, if so, from a STFT which also exhibits constant Heisenberg boxes in the time-frequency plane.

Thank you for helping in advance

A lifting scheme is a method for splitting a sequence of discrete samples into downsampled subsequences, so that you can predict a subsequence from the other, and update the former later, using possibly nonlinear predict and update operators. This corresponds to a generalization of one level of discrete wavelet transform, usually obtained fully linearly by convolution (first-generation wavelet).

Lifting wavelets were called second-generation wavelets, and proved to be a generalization of first-generation wavelets.

This being said, for question 1:

yes, a wavelet packet being iterations of several one level DWT on other subbands.

For question 2:

A lifting discrete wavelet packet can be  critically sampled, like the DWT, but also can be undecimated or oversampled if needed. And they are interesting questions about what kind of  stationary lifted wavelets can be built.

For question 3:

Generally, the Heisenberg boxes are considered to be of constant area. With the STFT, titles have the same shape, whereas with wavelet they are dilated/stretched in each direction. With a uniform wavelet packet decomposition (where the time-scale plot is regularly discretized), the time-scale bins can look similar to that of  the STFT. If the packet transform is critical, the result is often poorer with wavelets packets (because of the lack of redundnacy).

At the same redundancy, aliasing and constraints on wavelet shapes often yield less degrees of freedom that full-fledged oversampled filter banks, but the filter are generally easier to optimize. Nevertheless, if you want uniform frequency bin, wavelets are rarely very beneficial.


What is difference between the terms bit rate,baud rate and data rate? Let us consider any video signal. Bit rate of a video signal is the no of bits  per frame per second. Also, data rate is the no of bits per second. Baud rate is the number of signal changes per second.

But I really confused between these 3 terms. They look similar to each other.  So can anyone tell me whether these terms are same or different with the help any example?

Bit rate and data rate are somewhat ambiguous, and their exact definitions vary from one field or application to another.

Bit rate is often used to measure the information rate, since information is measured in bits. Sometimes it is used to measure the number of actual 1s and 0s transmitted per second. When the information is compressed (source-coded), the two definitions are equal.

Data rate can be used to mean the same as bit rate. Sometimes it is used to measure the rate at which data, as opposed to overhead, is transmitted. For instance, an audio compact disc system has a data rate of $2\cdot44100\cdot16=1411200$ bits (zeros and ones) per second, but the actual transmission rate is approximately double that, once you include framing, metadata, and FEC.

The baud rate (or symbol rate, or pulse rate) is the number of actual pulses (symbols) transmitted per second. If a pulse "carries" several bits, then the baud rate can be significantly smaller than the bit rate. This is important because the bandwidth required to transmit the information is determined by the baud rate.

So, consider an encoded video stream with 1000 bits per frame, at 24 frames per second, with 10% overhead per frame, and 4 bits per symbol:

From an information-theory perspective, the bit rate is $24\cdot 1000=24000$ (information) bits/sec.
You can also say, from a more practical perspective, that the bit rate required is $1.1\cdot 24 \cdot 1000 = 26400$ (ones and zeros) bits per second, but you could still define the data rate as $24000$ bits per second, since the remaining 2400 bits are overhead.
The pulse rate on the transmission link will be $26400/4=6,600$ pulses per second.


How do I calculate correlated color temperature correctly? How can I calculate the color temperature(CCT) of the light source illuminating a given raw image correctly?

The method I used is proposed by McCamy, and the algorithm was also described in the Wikipedia. The problem is that I don't know whether should I apply white balance process before calculating CCT.

If it should not, how do I get the conversion matrix that convert the current device-dependent RGB color space to CIE XYZ to calculate the algorithm correctly?

The problem is that I don't know whether should I apply white balance process before calculating CCT.

Definitely not. Afterwards, the light source would have a completely different color profile, depending on how you implement white balance, and hence, a different color temperature.

  If it should not, how do I get the conversion matrix that convert the current device-dependent RGB color space to CIE XYZ to calculate the algorithm correctly?

Get the color profile (which might be much more complex than a simple matrix) of your device; usually, that means getting the ICC file for your camera/scanner/sensor, and applying that to convert your image from the device color space to something canonical.

If your camera/scanner/sensor doesn't have such a profile, you will need to calibrate it first.


Pole/Zero Placement for Filter I'm attempting to design a FIR high pass filter than keeps signals above 200Hz and rejects signals below 60Hz with a sampling frequency of 500 samples/sec. This is my first time attempting this and I'm a little confused.

I started with looking at pole/zero placement. I'm not exactly sure how to go about this but this is what I did. I want a zero on the unit circle at the 60Hz point to minimize that frequency. To find this spot on the axis I did the following:

$$\frac{60}{\frac{500}{2}}\pi=0.24\pi$$

So I believe I need a zero at 0.24pi on the unit circle on the z-plane. Does this sound right? How would I then transform it in a difference equation in the sequence domain?

This is how you can try to design a short and simple FIR filter by hand. Note that this method is just supposed to be enlightening, a really useful filter can be designed using some software. If you want a zero at 60 Hz you indeed need to place it at an angle of $0.24\pi$ on the unit circle of the $z$-plane. The general relation between angle and frequency is

$$\phi=2\pi\frac{f}{f_s}$$

where $f$ is the desired frequency, and $f_s$ is the sampling frequency. Note that if the filter coefficients are real-valued, zeros (and poles) always occur in complex conjugate pairs. So if you have a zero at $z=z_0$ you must also have a zero at $z=z_0^*$. This just reflects the fact that the spectrum of a real-valued filter is symmetric: $H(\omega)=H^*(-\omega)$. Since you want to design a causal FIR filter, all poles are at the origin of the $z$-plane, and you only have control over the zeros.

So let's design an FIR filter of order $3$ (i.e. length $4$) with a zero at $z_0=e^{j0.24\pi}$, and - because it's a high pass filter - a zero at $z=0$. The zero at $z=0$ contributes a factor

$$H_1(z)=1-z^{-1}$$

to the total transfer function, and the zero at $z_0=e^{j0.24\pi}$ (and its complex conjugate counterpart at $z_0^*=e^{-j0.24\pi}$) contributes a factor

$$H_2(z)=(1-z_0z^{-1})(1-z_0^*z^{-1})=1-(z_0+z_0^*)z^{-1}+|z_0|^2z^{-2}=\\
=1-2\Re\{z_0\}z^{-1}+z^{-2}=1-2\cos(0.24\pi)z^{-1}+z^{-2}$$

because $z_0+z_0^*=2\Re\{z_0\}=2\cos(0.24\pi)$, and $|z_0|=1$ (the zero is on the unit circle). The total transfer function is

$$H(z)=H_1(z)H_2(z)=(1-z^{-1})(1-2\cos(0.24\pi)z^{-1}+z^{-2})=\\
=1-[1+\cos(0.24\pi)]z^{-1}+[1+\cos(0.24\pi)]z^{-2}-z^{-3}\tag{1}$$

From (1) you can read the filter coefficients:

$$h_0=1,\quad h_1=-[1+\cos(0.24\pi)]=-2.4579,\quad h_2=-h_1,\quad h_3=-h_0\tag{2}$$

Note that we just imposed the locations of the zeros but we didn't apply any normalization. A useful normalization for a highpass filter is to require $H(-1)=1$, i.e. the response at the Nyquist frequency (half the sampling frequency) is $1$. Since

$$H(-1)=\sum_n(-1)^nh[n]=h[0]-h[1]+h[2]-h[3]=2h[0]-2h[1]$$

we can define new normalized filter coefficients by

$$\hat{h}[n]=\frac{h[n]}{2(h[0]-h[1])},\quad n=0,1,2,3$$

The magnitude of the resulting frequency response of this filter looks like this (the phase is linear due to the symmetry of the coefficients):

You see that the three constraints (zeros at DC and at 60 Hz, and unity response at Nyquist) are satisfied. The transition from stopband to passband is of course not very steep but that's as good as it gets with an FIR filter of order 3.


How to Tell How Likely a Signal Is Present in Another One (Variance Unknown)? I know this is probably a simple question, but I haven't been able to find a satisfactory answer anywhere...

Say you have a time series signal of finite length N.  Call it $y[n]$.  It looks like a sine-gaussian perhaps but with some random effects.  The mean of zero, and no trend is present.

Now, say you receive this signal $y[n]$, and you have a set of K other 'candidate' signals, $x_{1}[n], x_{2}[n], ...,x_{K}[n]$.  You want to see what the most likely candidate is.  Is there a way to do this, and attach probabilities?  For example, candidate 1 has a 20% probability of being present, candidate 2 15%, etc. such that the percents add to 100%.

Some notes:

I want amplitude to matter.  If the candidate signal has a much smaller amplitude than $y[n]$, it should be less likely to be present, than another candidate which is otherwise exactly the same, but with a higher amplitude.
The variance of each data point in the signals (either $y[n]$ or $x[n]$) is unknown.  All we are given is what is mentioned above.  I know that in order to do something like a $\chi^{2}$ goodness of fit test (which has been suggested to me), something must be known about the variance of each data point in the signals.
The closest thing I have found is matched filtering, but how do I compute probabilities like I mentioned above?  Or is computing probabilities like that sort of the wrong answer to the question?
Coherence is related, but its more about how the signals change over time (from my limited understanding).  All the signals mentioned have a finite length N, and the signals are already matched in time (we only care about how similar they are at a particular time instant).  Time delays between them are irrelevant.

Thanks!! Any thoughts any of you have on this would be greatly appreciated!

Though the Matched Filter is the best tool detection of a known signals under AWGN it should work well here as well.

To say something about the probabilities the question is, do you know something about the energy of the received signals?
If you do, you should easily say something about the probabilities.

Pay attention that if the assumption is a signal could be within a different signal you need to apply a generalized Matched Filter (Multi Decision).


Build an inverse model for a train of gaussian pulses I have a stationary signal from a train of Gaussian pulses.
My sampling window is too wide (cannot be reduced). In the example 1 ms. So it is not possible to clearly distinguish one pulse from another.

The advantage is that I can resample the signal several times by moving sampling start with a delay. For example, I sample the signal 5 times, moving 0.2 ms in each one.

I need to build an inverse model that given the sampled signals (the 5 signals below), original signal can be reconstructed.

Is there any algorithm that solves this inverse problem?

The reconstruction can be done by commutating through each output resulting is a 5 kSps representation of the signal after having been passed through a 5 sample moving-average filter (5 tap FIR with unity gain coefficients). Thus any frequency content above 2.5 KHz will be aliased into the 0 to 2.5 KHz first Nyquist zone, and the overall result will be low pass filtered. Such a moving average filter would have its first null at 1 KHz. By reducing the time delay (and creating proportionally more outputs to commutate) the equivalent sampling rate of the resulting signal can be increased to reduce the aliasing but the first null of the resulting frequency response will always be at 1 KHz (the frequency response will approach a Sinc function as the time delay is reduced).

This is clear by observing the equivalent block diagram of the process below:

On each cycle of the 5 kSps clock, the input is fed into each 5 sample moving average filter ($MAF5$) with each filter input delayed by one sample. Once each filter has received 5 samples, the output is decimated by 5, thus providing a new sum over 5 samples, every 5 samples. Each decimated output is offset by one sample. Thus it is clear that each of these outputs is simply the output of the first $MAF5$ after each sample of the 5 kSps clock, as shown in the diagram below.

The magnitude of the frequency response for this process is shown in the figure below:

The process will result in distortion from both the aliasing from frequency content above the 2.5 KHz Nyquist Frequency and the frequency response of the 5 sample moving average filter. As stated in the introduction, reducing the delay and increasing the number of outputs will increase the effective sampling rate, thus reducing aliasing distortion, but the first null of the frequency response of the equivalent moving average filter will remain unchanged.


Doing resampling, how can we choose lower cutoff frequency for LPF filter? I would like to resample an audio signal.

Sampling frequency of an original audio signal is 50k Hz.

What I want to down sampling frequency is 30k Hz.

In this case, I followed the procedures as below using LabVIEW.

  '1. IIR low pass filter

  Butterworth / 6th / Lowpass / up to 15k Hz

  '2. Resample to 30kHz

Why I used LPF is to prohibit the anti-aliasing.

 from now on, I have some questions.

1. Do we need to do LPF for prohibiting anti-aliasing?

2. If we have to do, how can we choose the lower cutoff frequency?
   In my case, I choose as 15k Hz because this is the new Nyquist frequency.

Thank you in advance.

OR is there the standard procedure for resampling?

One thing to consider is the existing spectral content of the signal being downsampled.   If there is no significant energy above 15kHz, then the anti-alias filtering requirements are certainly relaxed since it won't have to do much for a 30kHz output sampling rate.

So where to put the filter cutoff frequency also depends on assumptions about the spectral content of the input signal as well as the transition band rolloff characteristics of the filter.  For a 30kHz sampling rate the 3dB cutoff point of the filter may need to be much lower than 15kHz to prevent aliasing from 15.5kHz or 16kHz (for example).

Much of this depends on your tolerance for the level of aliasing and at what frequencies.


What is anti-alias pre-filter for preventing aliasing after under-sampling? We know that the under-sampling results in aliasing and frequencies higher than half of the Nyquist rate is not distinguishable. I've a base band signal that I want to use the higher frequencies which are higher than the half of the Nyquist rate (Nyquist frequency) as well as low frequencies(all parts). I've a special process with this path:

$$\textrm{Input}{\longrightarrow}\boxed{\textrm{anti-aliasing pre-filter}}{\longrightarrow}\boxed{\textrm{decimate}}{\longrightarrow}\boxed{\textrm{FFT}}{\longrightarrow}\boxed{\textrm{tune on special part}\\{\textrm{of the signal}}}$$

The low-pass post-filter that people usually use as anti-aliasing filter removes the high frequencies that are of interest for me.
What is the digital or analog anti-aliasing pre-filter that I don't lose high frequencies.

I think you're looking for a free lunch that does not exist. Your original question and response to Peter K's answer suggests that you want to sample a signal that has both lowpass and highpass content, with the highpass content extending beyond the Nyquist frequency associated with your target sample rate. That is probably not going to work.

Given a sample rate $f_s$ (and real samples), you can only unambiguously represent frequencies on the interval $[0, \frac{f_s}{2})$. More generally, you can only represent a swath of bandwidth that is up to $\frac{f_s}{2}$ wide. Frequencies above the Nyquist rate alias down such that they appear to be located in this interval after you sample them. If you have a signal of interest that meets this bandwidth constraint, then you can use bandpass sampling techniques; basically you select a sample rate taking into account the center frequency and bandwidth of the desired signal. You allow the signal to alias in a "controlled" manner, so that it appears to be present in a contiguous portion of $[0, \frac{f_s}{2})$ after you sample it (perhaps with the spectrum inverted, but that is easily fixed).

This does not seem to line up well with what you seem to want. Your question indicated that you have lowpass content (i.e. content near zero frequency) that you want to preserve in addition to highpass content above $\frac{f_s}{2}$. In many cases, this is not going to be achievable without the highpass content aliasing down on top of the lowpass signal of interest after you sample. However, under certain conditions, you might be able to make this work. If:

The lowpass and highpass components are separated in frequency (i.e. there is a gap between the two regions where you don't care about preserving the signal's content),
You know the center frequency and bandwidth of the highpass portion (so it is more accurately termed "bandpass" instead),
And you have control over the sample rate,

Then you may be able to make it work. In that relatively special case, you would simply apply the bandpass sampling approach described before, except the sample rate must be selected with caution so that the higher-frequency content does not alias down into the portion of the band that the lowpass signal occupies.

Whether you would actually want to do this in a practical system is still an open problem. It's not clear specifically what you're trying to accomplish, or what the constraints in your application are. An alternate approach would be to separate the two signal components using analog filters (lowpass for one channel, highpass/bandpass for the other), then sample them independently. This could allow you to use a lower sample rate, commensurate with the bandwidths of each component.


Find $v_3$ such that $\{v_0, v_1, v_2, v_3\}$ is an orthogonal basis of $\mathbb{R}^4$? In this question, we consider the Hilbert space of vectors in $\mathbb{R}^4$ with one of its basis $\{v_0, v_1, v_2, v_3\}$.
Given the first three basis vectors
$$v_0=\left[\begin{array}{c}
1/2\\
1/2\\
1/2\\
1/2
\end{array}\right]v_1=\left[\begin{array}{c}
1/2\\
1/2\\
-1/2\\
-1/2
\end{array}\right]
v_2=\left[\begin{array}{c}
1/2\\
-1/2\\
1/2\\
-1/2
\end{array}\right]$$

How many possibilit(y/ies) are there for $v_3$, such that $\{v_0, v_1, v_2, v_3\}$ is an orthogonal basis of $\mathbb{R}^4$? Give one example of $v_3$ please.

let $v_3 = \left[\begin{array}{vector}
a\\
b\\
c\\
d\end{array}\right]$, then $<v_0, v_3>=0, <v_1, v_3>=0, <v_2, v_3>=0, <v_3, v_3>=1$,
reduce to

a+b=0
a+c=0
a-d=0

and $a^2+b^2+c^2+d^2 = 1$

so there are many solutions, one example can be

$ \left[\begin{array}{arr}
 \frac{1}{2}\\
 -\frac{1}{2}\\
 -\frac{1}{2}\\
 \frac{1}{2}\end
{array} \right]$


Suitable metrics for summarizing or visualizing the spatial activation of components after temporal Independent Component Analysis (ICA) I will first describe the Independent Component Analysis (ICA) steps, so that the question becomes clearer, and more relevant to similar questions.

Assume we have $N$ sensors distributed in space. From each sensor we obtain a time series of length $T$. We collect these time series as rows in an $N \times T$ matrix $X$. After performing ICA we obtain a decomposition of the form $X = DS$, where $D$ is an $N \times N$ matrix of weights, and $S$ is an $N \times T$ matrix whose rows are the independent sources. To find the contribution of the $i$-th component to the observed signal we then multiply the $i$-th column of the matrix $D$ with the $i$-th row of $S$. The result is an $N \times T$ matrix which we will denote as $X^{i} = D_{*,i}S_{i,*}$. Then we can write: $X = X^{1} + ... + X^{i} + ... + X^{N}$.

In many domains I see single-picture visualizations of the spatial activation of the $i$-th component from ICA. However, it is not clear to me what metric to extract from each row of $X^{i}$ or from matrix $D$ in order to do this. Is it a summary statistic of the rows of $X^{i}$, e.g. the energy of the signals (rows) of $X^{i}$, the weighting coefficients from the $i$-th column of $D$, or something different altogether?

Thank you.

As we know from the instantaneous mixture model of ICA, a signal mixture $X_{N,T}$ is composed of an unknown mixing matrix $D_{N,N}$, and an unknown signal matrix $S_{N,T}$, where:

$$
X = DS
$$

In ICA, we solve for $D_{N,N}$ and $S_{N,T}$ via maximization of absolute kurtosis. (As one method).

If we now look at the $i$th column of $D$, each element therein will now give us weight by which the $i$th signal gets weighted by in all $N$ mixtures. That is, each element $D_{n,i}$, (row $n$, column $i$), gives us signal $i$'s contribution to mixture $n$.

Those are called 'spatial weights' simply because in many ICA applications, we are de-mixing signals that sample a spatial field at different locations. Thus, each element of a column $D$ corresponds to a weight in space, that weighs signal source $i$.

In this way, we can show those weights in space as the heatmap above shows, and find the weight each signal gets weighted by, for a particular mixture, at different locations in space.


QAM for analog and digital waves I tried the following code for modulation of two analog waves. Can anyone point to how to modulate and demodulate for an analog and a digital wave? I tried by replacing the second signal (x2) with a digital pulse but I am not able to recover the analog wave after demodulation.

T = 1;
Fs = 44100;
t = 0:1/Fs:T;
Fc = 17000;

f1 = 100;
x1 = sin(2*pi*f1*t);
f2 = 200;
x2 = sin(2*pi*f2*t);

Cs = sin(2*pi*Fc*t);
Cc = cos(2*pi*Fc*t);

%%Modulation
mod_wav = x1.*Cs + x2.*Cc;

%%Demodulation
demod_x1 = mod_wav.*Cs;
demod_x2 = mod_wav.*-Cc;

%%filter
fpass = 10; fstop = 500; N = 30;
b = firls(N,[0 fpass fstop Fs/2]/Fs/2, [1 1 0 0], [1 1]);
filt_x1 = filter(b,1,demod_x1);
filt_x2 = filter(b,1,demod_x2);

In the above code, if i use butter filter, then it is properly working but not if i use FIR like above. Why is that? In a broader sense, when to use butter and and when to use FIR filter?

When I ran your code in Matlab I got the following warning message:

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND =  6.912349e-21.
> In firls at 175

It turns out that the problem is the "/Fs/2" part.  You are trying to normalize the frequencies from $0$ to $\frac{F_s}{2}$ into the range $0$ to $1$.  The problem is that Matlab divides the numbers in order, so even though what you want is /(Fs/2) what you are getting is /(Fs*2).  Changing the line to the following-

b = firls(N,[0 fpass fstop Fs/2]/(Fs/2), [1 1 0 0], [1 1]);

fixes the problem.


DFT of a series of RC exponentials Context: I'm trying use matlab to apply a single-pole filter to a time-domain ramp waveform that is generated by a sequence of time-shifted "RC steps" that are added together.

The time domain voltage waveform is

$V(t) = \sum_{k=0}^{N-1} V_{step}(1-e^{-(t-k\Delta t_{step})/\tau})u(t-k\Delta t_{step})$

where $V_{step}$ and $\Delta t_{step}$ are constants

Using the fft of a one-sided exponential decay, the unit-step, and the time-shift property of the fft I get a frequency domain representation of:

$\mathfrak{F}(V(t)) = V_{step}\left(\sum_{k=0}^{N-1}e^{-j\omega k\Delta t_{step}}\right)\left(\frac{1}{j\omega} +\pi\delta(\omega) - \frac{1}{\frac{1}{\tau}+j\omega}\right)$

So now I want to generate this complex-valued fft manually in matlab, multiply by a filter response, and inverse fft. To start, I'm checking the frequency response of the input and its ifft to make sure it looks right:

N=32;
Vstep=25e-3;
tstep=10e-12;
tau=5e-12;
ts=0.1e-12;

Nfft = 2^nextpow2(max([N*tstep/ts N*tau*5/ts])); %Get enough points for the whole ramp with at least 5 tau's per exponenetial
w=(0:Nfft-1)*2*pi/Nff;t
f=w/2/pi/ts;

timeshifts=sum(exp(-1i*(0:N-1)'*w*tstep));
step=[pi 1./(1i*w(2:end))];
expdecay=1./(1/tau+1i*w);

Vf=Vstep*timeshifts.*(step-expdecay);
vt=ifft(V);

So the frequency domain Vf looks reasonable in amplitude. However, when I take the inverse and plot vt, it is definitely not correct.

Where am I going wrong? I suspect there's something I'm missing with the fact that this is a DFT not fourier transform. Also I know there's some subtlety to the FT/DFT of the heaviside, particularly at w=0, and I know that the fft at 0 will just be the average (times N) and I'm not sure my script accomplishes this.

I do know that I can just start with the time domain and then fft it, but I'm rather curious now.

okay, for step 1, let's get the continuous time math down right and simplified to as much as possible.  i am changing some notational convention slightly to be more consistent with Electrical Engineering practice:

$$v(t) = \sum_{k=0}^{N-1} V_\mathrm{step}(1-e^{-(t-k\Delta t_\mathrm{step})/\tau})u(t-k\Delta t_\mathrm{step})$$

$$V(j\omega) =\mathscr{F}\Big\{v(t)\Big\} = V_\mathrm{step}\left(\sum_{k=0}^{N-1}e^{-j\omega k\Delta t_\mathrm{step}}\right)\left(\frac{1}{j\omega} +\pi\delta(\omega) - \frac{1}{\frac{1}{\tau}+j\omega}\right)$$

there is a well known closed form for the summation of a geometric series:

$$ \sum_{k=0}^{N-1} x^k = \frac{x^N -1}{x-1} $$

so

$$\begin{align}
 \sum_{k=0}^{N-1} e^{-j\omega k\Delta t_\mathrm{step}} &= \sum_{k=0}^{N-1} (e^{-j\omega \Delta t_\mathrm{step}})^k \\ \\
 &= \frac{(e^{-j\omega \Delta t_\mathrm{step}})^N -1}{e^{-j\omega \Delta t_\mathrm{step}}-1} \\ \\
 &= \frac{e^{-j\omega N \Delta t_\mathrm{step}}-1}{e^{-j\omega \Delta t_\mathrm{step}}-1} \\ \\
 &= \frac{e^{-j (\omega/2) N \Delta t_\mathrm{step}}(e^{-j (\omega/2) N \Delta t_\mathrm{step}}-e^{j (\omega/2) N \Delta t_\mathrm{step}})}{e^{-j(\omega/2) \Delta t_\mathrm{step}}(e^{-j(\omega/2) \Delta t_\mathrm{step}}-e^{j(\omega/2) \Delta t_\mathrm{step}})} \\ \\
 &= e^{-j (\omega/2) (N-1) \Delta t_\mathrm{step}}\frac{(e^{j (\omega/2) N \Delta t_\mathrm{step}}-e^{-j (\omega/2) N \Delta t_\mathrm{step}})/(2j)}{(e^{j(\omega/2) \Delta t_\mathrm{step}}-e^{-j(\omega/2) \Delta t_\mathrm{step}})/(2j)} \\ \\
 &= e^{-j (\omega/2) (N-1) \Delta t_\mathrm{step}}\frac{\sin(\omega N \Delta t_\mathrm{step}/2)}{\sin(\omega \Delta t_\mathrm{step}/2)} \\ \\
\end{align}$$

so i think your final expression in continuous frequency of the spectrum should be

$$V(j\omega) = V_\mathrm{step}\frac{\sin(\omega N \Delta t_\mathrm{step}/2)}{\sin(\omega \Delta t_\mathrm{step}/2)} \left(\frac{1}{j\omega} +\pi\delta(\omega) - \frac{\tau}{1+j\omega \tau}\right) e^{-j (\omega/2)(N-1)\Delta t_\mathrm{step}} $$

does that look right to you?


Extended Kalman Filter for Linear Systems with Non Linear Measurements I'm successfully using an Extended Kalman Filter for object tracking. My state vector ($x, y, v_x, v_y$) needs to be in cartesian coordinates. The measurement data is transmitted in polar coordinates.

The state transition equations (and thus the state transition matrix) are linear:
\begin{align}
x_{k+1} &= x_{k} + v_{x,k} \cdot \Delta t \\
y_{k+1} &= y_{k} + v_{y,k} \cdot \Delta t \\
v_{x,k+1} &= v_{x,k} \\
v_{y,k+1} &= v_{y,k}.
\end{align}

So the non-linear relation of the measurement vector in polar coordinates to my state vector in cartesian coordinates is the only reason why I used an EKF.

A simpler approach could include a "classic" transformation of my measurement data from polar to cartesian coordinates before feeding it into a standard Kalman Filter.

Assuming that I transform my measurement noise covariance matrix accordingly: Would this have an impact on the filter performance (e.g. because the measurement noise can't be assumed as gaussian any more)?

When you start comparing sub optimal approaches, it usually comes down to trial and error.

A KF essentially operates on means and variances. An EKF does the same but the nonlinear measurement and or states necessitate some additional considerations.

The first thing one needs to determine is if the transformation results in a unimodal probability distribution where there is a well defined mean and variance.  If not, a particle filter will be a better solution. If the unimodal assumption is realistic, then the bias associated with how you map or update a mean is the next problem. In some cases a linearization is good enough and sometimes the transformation needs to be used and that includes things like Runge Kutta ODE solvers, which most text books fail to mention. In your case, a measurement nonlinearity can either use a transformation or a linearization based on a Taylor expansion to propagate a mean.

The final problem is how to propagate the variance.  Most EKF will use a Taylor Expansion most typically, just the first derivative, but there are second and higher order filters as well.  These enter as complexity considerations as well. The text edited by Gelb has a good number of covariance update heuristics in the latter chapters, including those for discontinuous transformation.

Gelb, Arthur, ed. Applied optimal estimation. MIT press, 1974.

The UKF is an attractive approach to covariant update.

Step size is also important. Your measurements may not be at regular time steps.

In summary, sub optimal is less than optimal but which suboptimal is better is going to depend on the application.  Essentially you need to try some alternatives and judge which is better.  There are a lot of alternatives and which is best isn’t easily predictable.


How does adaptive Huffman coding work? Huffman coding is a widely used method of entropy coding used for data compression. It assumes that we have complete knowledge of a signal's statistics. However, there are versions of Huffman coding that are used with streaming media and cannot possibly know everything about the signal's statistics. How do these adaptive Huffman encoders work?

The Wikipedia article has a pretty good description of the adaptive Huffman coding process using one of the notable implementations, the Vitter algorithm. As you noted, a standard Huffman coder has access to the probability mass function of its input sequence, which it uses to construct efficient encodings for the most probable symbol values. In the prototypical example of file-based data compression, for example, this probability distribution can be calculated by histogramming the input sequence, counting the number of occurrences of each symbol value (symbols could be 1-byte sequences, for example). This histogram is used to generate a Huffman tree, like this one (taken from the Wikipedia article):

The tree is arranged by decreasing weight, or probability of occurrence in the input sequence; leaf nodes at the top represent the most probable symbols, which therefore receive the shortest representations in the compressed data stream. The tree is then saved along with the compressed data and is subsequently used by the decompressor later to regenerate the (uncompressed) input sequence again. As one of the early entropy code implementations, standard Huffman coding is fairly straightforward.

The adaptive Huffman coder's structure is quite similar; it uses a similar tree-based representation of the input sequence's statistics to select efficient encodings for each input symbol value. The main difference is that, as a streaming implementation of the algorithm, no a priori knowledge of the input's probability mass function is available; the sequence's statistics must be estimated on the fly. If one is to use the same Huffman encoding scheme, this means that the tree used to generate each symbol's encoding in the compressed stream must be built and maintained dynamically as the input stream is processed.

The Vitter algorithm is one way of accomplishing this; as each input symbol is processed, the tree is updated, maintaining its characteristic of decreasing probability of symbol occurrence as you move down the tree. The algorithm defines a set of rules for how the tree is updated over time, and how the resulting compressed data is encoded in the output stream. As the input sequence is consumed, the tree's structure should represent a more and more accurate description of the input's probability distribution. In contrast to the standard Huffman coding approach, the decompressor does not have a static tree to use for decoding; it must perform the same tree-maintenance functions continuously during the decompression process.

In summary: The adaptive Huffman coder operates very similarly to the standard algorithm; however, instead of a static measurement of the entire input sequence's statistics (the Huffman tree), a dynamic, cumulative (i.e. from the first symbol to the current symbol) estimate of the sequence's probability distribution is used to encode (and decode) each symbol. In contrast to the standard Huffman coding approach, the adaptive Huffman algorithm requires this statistical analysis at both the encoder and decoder.


How to build a median filter function? What algorithms are there for 1-d median filtering (sliding median, rolling median, running median), similar to MATLAB's medfilt1? Of interest would be a reference implementation written from scratch, preferentially for MATLAB.

here is code for a "rolling median" copied from the reference of the reference i made above.  i believe the computational cost of the alg (per sample) is $O(\log_2(N))$.

//Copyright (c) 2011 ashelly.myopenid.com under <http://w...content-available-to-author-only...e.org/licenses/mit-license>

#include <stdlib.h>

//Customize for your data Item type
typedef int Item;
#define ItemLess(a,b)  ((a)<(b))
#define ItemMean(a,b)  (((a)+(b))/2)

typedef struct Mediator_t
{
   Item* data;  //circular queue of values
   int*  pos;   //index into `heap` for each value
   int*  heap;  //max/median/min heap holding indexes into `data`.
   int   N;     //allocated size.
   int   idx;   //position in circular queue
   int   ct;    //count of items in queue
} Mediator;

/*--- Helper Functions ---*/

#define minCt(m) (((m)->ct-1)/2) //count of items in minheap
#define maxCt(m) (((m)->ct)/2)   //count of items in maxheap

//returns 1 if heap[i] < heap[j]
int mmless(Mediator* m, int i, int j)
{
   return ItemLess(m->data[m->heap[i]],m->data[m->heap[j]]);
}

//swaps items i&j in heap, maintains indexes
int mmexchange(Mediator* m, int i, int j)
{
   int t = m->heap[i];
   m->heap[i]=m->heap[j];
   m->heap[j]=t;
   m->pos[m->heap[i]]=i;
   m->pos[m->heap[j]]=j;
   return 1;
}

//swaps items i&j if i<j;  returns true if swapped
int mmCmpExch(Mediator* m, int i, int j)
{
   return (mmless(m,i,j) && mmexchange(m,i,j));
}

//maintains minheap property for all items below i/2.
void minSortDown(Mediator* m, int i)
{
   for (; i <= minCt(m); i*=2)
   {  if (i>1 && i < minCt(m) && mmless(m, i+1, i)) { ++i; }
      if (!mmCmpExch(m,i,i/2)) { break; }
   }
}

//maintains maxheap property for all items below i/2. (negative indexes)
void maxSortDown(Mediator* m, int i)
{
   for (; i >= -maxCt(m); i*=2)
   {  if (i<-1 && i > -maxCt(m) && mmless(m, i, i-1)) { --i; }
      if (!mmCmpExch(m,i/2,i)) { break; }
   }
}

//maintains minheap property for all items above i, including median
//returns true if median changed
int minSortUp(Mediator* m, int i)
{
   while (i>0 && mmCmpExch(m,i,i/2)) i/=2;
   return (i==0);
}

//maintains maxheap property for all items above i, including median
//returns true if median changed
int maxSortUp(Mediator* m, int i)
{
   while (i<0 && mmCmpExch(m,i/2,i))  i/=2;
   return (i==0);
}

/*--- Public Interface ---*/

//creates new Mediator: to calculate `nItems` running median.
//mallocs single block of memory, caller must free.
Mediator* MediatorNew(int nItems)
{
   int size = sizeof(Mediator)+nItems*(sizeof(Item)+sizeof(int)*2);
   Mediator* m=  malloc(size);
   m->data= (Item*)(m+1);
   m->pos = (int*) (m->data+nItems);
   m->heap = m->pos+nItems + (nItems/2); //points to middle of storage.
   m->N=nItems;
   m->ct = m->idx = 0;
   while (nItems--)  //set up initial heap fill pattern: median,max,min,max,...
   {  m->pos[nItems]= ((nItems+1)/2) * ((nItems&1)?-1:1);
      m->heap[m->pos[nItems]]=nItems;
   }
   return m;
}

//Inserts item, maintains median in O(lg nItems)
void MediatorInsert(Mediator* m, Item v)
{
   int isNew=(m->ct<m->N);
   int p = m->pos[m->idx];
   Item old = m->data[m->idx];
   m->data[m->idx]=v;
   m->idx = (m->idx+1) % m->N;
   m->ct+=isNew;
   if (p>0)         //new item is in minHeap
   {  if (!isNew && ItemLess(old,v)) { minSortDown(m,p*2);  }
      else if (minSortUp(m,p)) { maxSortDown(m,-1); }
   }
   else if (p<0)   //new item is in maxheap
   {  if (!isNew && ItemLess(v,old)) { maxSortDown(m,p*2); }
      else if (maxSortUp(m,p)) { minSortDown(m, 1); }
   }
   else            //new item is at median
   {  if (maxCt(m)) { maxSortDown(m,-1); }
      if (minCt(m)) { minSortDown(m, 1); }
   }
}

//returns median item (or average of 2 when item count is even)
Item MediatorMedian(Mediator* m)
{
   Item v= m->data[m->heap[0]];
   if ((m->ct&1)==0) { v= ItemMean(v,m->data[m->heap[-1]]); }
   return v;
}

/*--- Test Code ---*/
#include <stdio.h>
void PrintMaxHeap(Mediator* m)
{
   int i;
   if(maxCt(m))
      printf("Max: %3d",m->data[m->heap[-1]]);
   for (i=2;i<=maxCt(m);++i)
   {
      printf("|%3d ",m->data[m->heap[-i]]);
      if(++i<=maxCt(m)) printf("%3d",m->data[m->heap[-i]]);
   }
   printf("\n");
}
void PrintMinHeap(Mediator* m)
{
   int i;
   if(minCt(m))
      printf("Min: %3d",m->data[m->heap[1]]);
   for (i=2;i<=minCt(m);++i)
   {
      printf("|%3d ",m->data[m->heap[i]]);
      if(++i<=minCt(m)) printf("%3d",m->data[m->heap[i]]);
   }
   printf("\n");
}

void ShowTree(Mediator* m)
{
   PrintMaxHeap(m);
   printf("Mid: %3d\n",m->data[m->heap[0]]);
   PrintMinHeap(m);
   printf("\n");
}

int main(int argc, char* argv[])
{
   int i,v;
   Mediator* m = MediatorNew(15);

   for (i=0;i<30;i++)
   {
      v = rand()&127;
      printf("Inserting %3d \n",v);
      MediatorInsert(m,v);
      v=MediatorMedian(m);
      printf("Median = %3d.\n\n",v);
      ShowTree(m);
   }
}


Why is an interleaver often placed between an outer block code and an inner convolution code? It seems that when doing FEC with a concatenated-code, an interleaver is often placed between the outer block code and the inner convolution code.  The explanations I've seen say that this is because the convolution decoder often has errors in bursts.

However, if we assume that we're using a Reed-Solomon decoder with 8-bit symbols, then spreading the bursts onto multiple blocks would seem to increase errors rather than decrease them.  What am I missing that explains why this interleaving is good?

Rough idea of the decode and encode process:

Encode:

$$\boxed{\textrm{Data}}{\longrightarrow}\boxed{\textrm{RS Encoder}}{\longrightarrow}\boxed{\textrm{Interleaver}}{\longrightarrow}\boxed{\textrm{Convolution Encoder}}{\longrightarrow} \boxed{\textrm{Modulator}}$$

Decode:

$$\boxed{\textrm{Demodulator}}{\longrightarrow}\boxed{\textrm{Convolution Decoder}}{\longrightarrow}\boxed{\textrm{De-interleaver}}{\longrightarrow}\boxed{\textrm{RS Decoder}}{\longrightarrow}\boxed{\textrm{Data}}$$

As the other answers mention, the use of FEC results in post-decoding errors occurring
in bursts. Indeed, this happens regardless of whether the code is a convolutional code
or a block code.  With a $(n,k)$ block code, the decoder output ($k$ bits)
from the decoding of one
received word is (hopefully with high probability) completely correct, or it has
an unknown number of errors in it that can be regarded as a burst error of length $k$.
With a convolutional code, the decoder output is mostly correct as the decoder
finds the correct path through the trellis, but occasionally the decoder's
chosen path deviates from the correct path and later rejoins the correct
path during which time there is a burst of errors. In contrast to
block codes, there can be multiple
isolated burst errors (of variable lengths) in a single transmission
using a convolutional code. Also, the number of data bits in a single
transmission is far larger than the typical values of $k$ for a
block code.

The idea behind a concatenated
coding scheme is that with an inner code suited to the physical channel,
and an outer code over a very large symbol alphabet, we can make
the burst errors in the inner decoder output look like a single
symbol error to the outer code.  This is important because the
outer code should be a very high rate code because the net rate
is product of the inner code rate (more or less determined by
the channel and the link budget) and the outer code rate.
Unfortunately, outer codes over very large alphabets are very
expensive to implement, and so interleaved Reed-Solomon codes
over smaller alphabets (often GF$(2^8)$ ) are used (with
interleaving at the symbol
level, as Jim Clay points out). Because of the interleaving, the burst
errors in the inner decoder output become single symbol (byte) errors in the
received words of the interleaved Reed-Solomon code.

All the above is mostly a rehash of what the answers by Bryan and
Jim Clay have already said, but I wish to point out the following.

  Interleaved Reed-Solomon codewords can be decoded much more
  efficiently and with smaller delay if they are not de-interleaved
  first.

A Reed-Solomon decoder that can decode interleaved codewords is different from
the off-the-shelf standard Reed-Solomon
decoders that are available, and the use of such a decoder might not be
feasible if the development team does not have control of this aspect of
the design. But, if such a decoder is used, the de-interleaver
can be moved from the between the inner decoder and outer decoder to
just after the outer decoder. The de-interleaver is also smaller since
it has to deinterleave a $K\times L$ array instead of a $N \times L$
array (for a $(N.K)$ Reed-Solomon code interleaved to depth $L$.

  If a delay-scaled Reed-Solomon encoder is used along with the
  Reed-Solomon decoder for interleaved codewords described above, the interleaver
  at the transmitter and the de-interleaver at the receiver
  can be eliminated entirely.

The output of the delay-scaled
encoder (see also
this paper
which is unfortunately behind a paywall)
is a set of interleaved Reed-Solomon codewords, but
is not the same sequence of bytes that one would get from
doing a standard Reed-Solomon encoding followed by interleaving.
So, no further interleaving is necessary.
The output of the interleaved Reed-Solomon decoder is the
same byte sequence in the same order that went into the delay-scaled
Reed-Solomon encoder, and so no de-interleaving is necessary
at the decoder, either.


Detecting a three-tone signal in real-time I'm new to digital signal processing, and I've been looking into FFT, filtering, Matlab etc. What I'd like is to build a program which will detect when a certain BEEP is transmitted through a trunked voice channel (it signifies an impending critical message) being received via SDR. Through FFT I have ascertained it consists of three frequencies (500Hz, 600Hz and 700Hz). In essence -- somewhat like a phone call with DTMF -- I'd like to monitor audio for the presence of this beeeeep tone and have the software recognize when this is occurring.

My concern is that voice transmissions will cause these three frequencies to peak, making it difficult to determine whether it is due to the alert beep or somebody speaking.

I believe this might be accomplished through some means of measuring the length of the beeps, however this may be a misunderstanding through my inexperience in the field.

I'd appreciate any advice, and if my question isn't directly answerable, please push me in the right direction.

Edit: Here is a temporal frequency analysis of the series from WavePad

The approach I'm familiar with has been used for quite some time in broadcast situations. You mentioned your concern with voice signal "peaking" the frequencies of interest. From that statement, I am assuming that you can transmit your signal in-between periods of voice activity.

Here is the method:

On the generation side:

arrange for the signal to play for twice as long as the analysis.

On the analysis side:

Setup non-overlapped FFTs of a length that is 1/2 of the stimulus. This ensures that one of your FFTs will eventually contain an entire signal.
On each FFT, compute the power sum of the tone bins (including adjacent bins as a result of windowing), the power sum of the non-tone bins, and ratio the result to provide a SNR.
If the ratio is above some threshold, then you have detected your signal. You can determine your threshold using your prior knowledge of the noise floor and transmitted tone levels.

A couple of things you can do to improve:

Adjust the phases of the generated signal to reduce the crest factor.
Choose the 3 signal frequencies and/or the FFT length so that the frequencies fall in the centers of the FFT bins. If your generator and analyzer are synchronous, this would allow you to avoid windowing for FFT altogether.


Feature Extraction of Insect Sounds First of all, forgive me if I sound stupid. I am still learning Audio Processing.
So I am trying to make a machine learning model for detecting a specific insect through sound.

Here is a sample raw audio waveplot using librosa library:

After some preprocessing, applying High- pass filter and normalize it here is the output:

My questions are:

The next step would be to extract features from this, do I have to frame the audio manually and extract the features from that specific frame? Like this?

How would I remove the silence where the sound of the insects is not present?
Or do I have to analyze it directly for the whole duration of the sound? I have tested Librosa's Library in generating Spectrogram, here is the output:

Any good resources you can recommend for beginners?

Thank you so much.

Typically, for analyses like this, you do overlapping windows with a hop size that is fixed fraction of the window size.  So, for example, you might start by doing analyses of 1024 frames per window, and move the window forward 256 samples per analysis.  Depending on how the data looks then, you may or may not want to apply a window function (like Hamming or the like) to each window.

But since you're basically trying to remove noise, it might be simpler to take an example of the noise itself, and then try to reduce that noise from the original sample.  Take the FFT of the noise, subtract the frequency bins of the noise from the frequency bins of the sample, and run an inverse FFT to produce the final filtered result.  If you take this approach, make sure to only do FFTs at half the sample rate.


Computing the center of mass of a signal using the Fourier transform I am having difficulty with a homework problem which asks:

I have some ideas in mind but I have no clue as to whether they are correct or not. Below is my attempt:

The first thing you do is correct. You can calculate the sum of all the values of a signal evaluating its DTFT at $\omega=0$.

When you try to express the numerator with the DTFT though, you make a pretty clever thing but it is not correct, because you are evaluating the DTFT at $\omega=\ln\left( n^{\frac{1}{j2\pi n}}\right)$, which is not always real-valued. Remember that $\omega\in\mathbb{R}$.

To find how the DTFT of $x(n)$ and the numerator of $c$ are related, you can use the following property:

$$x(n)\xrightarrow{\mathscr{F}}X(\omega) \implies nx(n)\xrightarrow{\mathscr{F}}j\frac{dX(\omega)}{d\omega}$$

Using this property and your approach for the denominator (evaluating the transform at $\omega=0$), you should get the desired result.


Analytically designing a notch-filter for specified frequency 50 Hz I tried to analytically construct a causal notch filter that stops frequencies at 50Hz and I thought if zeros are more than poles it will be non-causal so I chose equal amount of poles and zeros resulting in

$$\dfrac{\left(z-e^{100\pi i} \right)\left(z-e^{-100\pi i} \right)}{z^2} $$

Is this correctly designed and causal?

You've designed a causal filter with a notch at $\omega_0=100\pi$. But the result is probably not what you want. Note that you've designed an FIR (finite impulse response) filter. Its frequency response has a large overshoot towards high frequencies.

What you actually want is an IIR filter, with poles away from the origin of the complex $z$-plane. A simple way to design a notch filter is to place the zeros at the desired notch frequency (as you did), and place poles at the same angle but slightly inside the unit circle. If $0<r<1$ is the pole radius, your transfer function looks like this:

$$H(z)=\frac{(z-e^{j\omega_0})(z-e^{-j\omega_0})}{(z-re^{j\omega_0})(z-re^{-j\omega_0})}=\frac{z^2-2z\cos(\omega_0)+1}{z^2-2zr\cos(\omega_0)+r^2}\tag{1}$$

You can apply scaling to the transfer function $(1)$ such that its gain at DC equals 1.

Finally, note that you have to take the sampling frequency into account. The notch frequency $\omega_0$ in $(1)$ is normalized. If $f_0$ is the desired notch frequency in Hertz, and if $f_s$ is the sampling frequency in Hertz, then

$$\omega_0=2\pi\frac{f_0}{f_s}$$

As an example, take $f_0=50\text{ Hz}$, $f_s=1000\text{ Hz}$, and $r=0.98$. The figure below shows the transfer function of the FIR notch filter you proposed (top), and the IIR transfer function according to Equation $(1)$ (bottom):


MFCC in Speech Recognition and ANN I have a dataset with over 2000 examples of digits (0 to 9) recorded on wav files. I'm trying to process these files with MFCC so that later I could train an Artificial Neural Network (ANN).

My problem is that ANNs requires, each sample of data, to have the same size, and every single tool I tried (HTK, Matlab and others..), generates outputs of different sizes.

This is probably due to the fact that some digits take more time to be pronounced than others. But how could I deal with this? Should I crop then in order to force them to have the same size? Or should I concatenate all of them with '0's, which would also make them with the same size?

I got stuck with the implementation, any suggestion is welcome

You might try padding the shorter (and even the longer) speech samples with a wide range of different (but realistic) background noise environments to create a larger training set (all to an equal training input size).  This might allow the DNN to be trained to be more robust at inference in noise.

You might even want to add variations of noise and a range of slight time-pitch variations to all your dataset samples to enlarge your training set.