DSP filter. Is prewarping performed when discretizing using Forward/Backward? I am trying to derive the coefficients used for a IIR implementation for the lowpass portion of a SVF filter. I've seen finished expressions for the coefficients (smith (p.8) and victor), their discretizized transfer function looks something like (not super important for my question)

$$
H(z) = \frac{k_f^2 \, z^{-1}}{1-(2-k_f(k_f+k_q))z^{-1} + (1-k_f k_q)z^{-2}} \, ,
$$
with $k_q = \frac{1}{Q}$ and (next is central to my question) $k_f = 2 \, \sin (\pi f_c T) $ .

I would like to derive these parameters myself, which I have attempted, but I am having trouble with the parameters involving the cutoff frequency $f_c$!

This is my starting point, (from analysing SVF circuit), I get the following transfer function for the lowpass portion,

$$
H(s) = \frac{\Omega_0^2}{s^2 + \frac{\Omega_0}{Q}s + \Omega_0^2} \, .
$$

Next step is performing the discretization. I did this using the bilinear transform ($z = e^{sT}$), using from what I gather the Tustin/trapezoidal approximation ($s = \frac{1}{T} \log(z) \approx \frac{2}{T} \frac{z-1}{z+1})$.

Prewarping

This transform requires prewarping and was done as follows,

$$ \Omega_0 = \tfrac{2}{T} \tan \left( \tfrac{\omega_0}{2} \right) \,, $$

where $\Omega_0 = 2 \pi f_c$ and $f_c$ is my wanted cutoff frequency.

Now, already I can see that my derivation is going to result in a different expression than what the two links above are getting. Their frequency parameters shows up in expressions that look like $ f = 2 \sin(\pi f_c T) $

My Question

Basically, where is the $\sin()$ term coming from?

Could it be that they use Euler forward (or backward) as an approximation?

$$ s = \frac{1}{T} \log(z) \approx \frac{z - 1}{T} $$

And then done some prewarping for that case? I have looked around the web and can not find any mention of prewarping done for euler backward/froward. Is that usually done?

That would give me
$$ \begin{split}
s = j \Omega &= \frac{z-1}{T} \\
j \Omega &= \frac{e^{j \omega}-1}{T} \\
j \Omega &= \frac{e^{j \omega/2}-e^{-j \omega/2}}{T e^{-j \omega/2}} \\
j \Omega &= \frac{j 2 \sin(\omega/2)}{T e^{-j \omega /2}} \\
 \Omega &= \frac{2 \sin(\omega/2)}{T e^{-j \omega /2}} \\
\end{split} $$

Not sure what I am missing.

--

Thank you very much for any help regarding this!

(After some helpful comments, I think I see it a bit more clearly. Here is an attempt at answering my own question, please let me know if I am completely off track.)

Answer:

In my question I must have misunderstood from where the top expression $H(z)$ comes from.
I thought that such an expression was the result of a direct derivation from an analog circuit using the bilinear transform.

As Robert comments to the original question, $H(z)$ as I've written it, is the transfer function for the Hal Chamberlin SVF. And is not the result of a direct derivation from an analog transfer function using the bilinear transform (also pointed out by Matt).

Also, to answer my other question as far as I know there is no pre-warping done for Euler forward/backward and the term $\sin(\pi f_c T)$ does not come up because of this.

Where does the $\sin(\pi f_c T)$ term come from?

Bluntly, I think from an approximation. An approximation that is baked into to what we call the Hal Chamberlin SVF (presented in "Musical Applications of Microprocessors").

By using Roberts fantastic EQ cookbook, which has transfer functions and expressions for the coefficients, but which are derived from analog prototypes, I can compare and see where the approximation is being made.

(I received a nice guide that walks through the Hal Chamberlin SVF its transfer function and how it compares to transfer functions derived using the Bilinear transform and later the Impulse Invariant mapping. Link here: "Exact Bilinear Transform or Impulse Invariant mappings of Q and resonant frequency from analog prototype to Hal Chamberlinʼs State Variable Filter (SVF)" by Robert Bristow-Johnson. I'll also paraphrase the part which explains the $\sin{()}$ term.)

Starting from (note: not derived via bilinear transform)

$$
H(z) = \frac{k_f^2 \, z^{-1}}{1-(2-k_f(k_f+k_q))z^{-1} + (1-k_f k_q)z^{-2}} \, ,
$$

I can compare the coefficients in the denominator with the coefficients that were derived from an analog lowpass biquad which are presented in the EQ cookbook. From this we see that

$$
  \frac{a_2}{a_0} = \frac{1 - \frac{\sin(\omega_c)}{2 Q}}{1+\frac{\sin(\omega_c)}{2 Q}} = 1 - k_f k_q
$$
$$
  \frac{a_1}{a_0} = \frac{ - 2 \cos(\omega_c)}{1+\frac{\sin(\omega_c)}{2 Q}} = k_f^2 + k_f k_q - 2
$$

where $\omega_c = 2 \pi f_c T$.

Using these two and solving for $k_f$ gives us,

$$
k_f = \frac{2 \sin(\frac{\omega_c}{2})}{\sqrt{1+\frac{\sin(\omega_c)}{2 Q}}} \, ,
$$

If we assume $Q > \frac{1}{2}$ and that $\omega_c \ll \pi$ (cutoff frequency is small compared to Nyquist) we can approximate,

$$
k_f \approx 2 \sin\left(\frac{\omega_c}{2} \right) = 2 \sin(\pi f_c T) \, ,
$$

and we have the $\sin(\cdot)$ term I was looking for.


MATLAB: How to filter a signal using Thiran Fractional Delay Filter I have generated a Thrian fractional delay filter using (thiran function) in MATLAB

h=thiran(0.3,1)

where 0.3 is the time delay and 1 is the sample time. The output is a discrete time transfer function, the generated filter is all-pass, infinite impulse response (IIR) filter.

The problem is that I could not apply it to a signal x(t). So I would like to know how to filter a signal with the generated thiran filter in MATLAB, what should I use to filter a signal with thiran filter? I can generate thiran filter using (thiran) function and the result will be a discrete time tf, but how to use it to filter a signal. I have tried to use filter, conv and tf functions but they all did not work:

filter and conv functions did not work because they require polynomial coefficients or numerator and denominator coefficients b and a, whereas thiran function generates a discrete transfer function.

[num,den] = tf(d)

function did not work because it can only convert a filter transfer function which is configured by designfilt function, so in this case d cannot be a thiran filter, d can only be configured by designfilt function

try this:

h=thiran(0.3,1);
num=  cell2mat(h.Numerator);
denum=cell2mat(h.Denominator);
output=filter(num,denum,x)

Hope that helps


How to detect the maximum resolvable spatial frequency of camera? I am trying to calculate the minimum line pixel width that can be distinguished from noise as shown in the camera test chart in Figure 1 where the thinner lines on the left are getting more and more blurry and mixed with noise thus become unresolvable. I took the picture of the camera quality test chart using an LG V30 smartphone. The ideal image should look like Figure 2.

Does anyone have an idea about how to mathematically calculate the frequency or pixel size at which the lines are unresolvable or where lines becomes visible as noise?

I figured that the Fourier transform can be used (which is the bottom plot of Figure 1, but still struggling to know how).

Figure 1

Figure 2

The image that you have obtained is already capturing the Modulation Transfer Function (MTF) of the optical system under test. It does not need to be further passed via an extra Discrete Fourier Transform (DFT) step. The spatial resolution limit is a matter of overall specification.

The test pattern's spatial resolution is measured in lines per millimeter and it varies in (a typical) $x$ direction (the width of the image) increasing towards the right.

As the resolution is increasing, the contrast (between the dark and light areas of the pattern) is decreasing.

Therefore, the simplest way to obtain the MTF is to create a scatter plot whose $x$ axis is spatial resolution in 1, 2,3,... lines/mm and the $y$ axis is the contrast (or modulation factor for that particular frequency) as $\frac{I_{max}-I_{min}}{I_{max}+I_{min}}$.

BUT, as you can see, the difference between $I_{max}, I_{min}$ in the test image is dropping to zero (blue trace in your diagram). Therefore, the test pattern (image) itself does not represent a $1$ input to the optical system across all spatial frequencies.

Therefore, it would be best to calculate the MTF of the test pattern and then the MTF of the system under test and present both of them on the same diagram if you want to depict how close to the "ideal" the optical system under test really is. This is the point of the MTF.

To find the smallest object that this optical system can resolve, you need to translate the size of that object into lines/mm and check that the MTF at that limit is above the noise floor of the sensing element.

Hope this helps.


Detection of a Circle in Noisy Image Data I have an image that looks like the one below:

I'm trying to find the radius (or diameter) of the circle.
I have tried using circular Hough transform (via matlab's imfindcircles(bw,[rmin rmax],'ObjectPolarity','bright')) , and by fitting to a circle or an ellipse (home made function that works pretty well for less noisy data, see below).

I've also tried some image processing to get a clearer circle, for example, see below:

se = strel('disk', 2);
bw = imdilate(bw, se);
bw = bwareaopen(bw,100000);
bw =  edge(bw);

However, when I feed the processed image to either techniques (Hough and circle\ellipse fitting) neither of them manage to detect the circle in a decent manner.

Here's a code snippet of the circle finder I wrote (matlab)
       [row col]=find(bw);
       contour = bwtraceboundary(bw, row(1), col(1)], 'N', connectivity, num_points);

    x = contour(:,2);
    y = contour(:,1);

    % solve for parameters a, b, and c in the least-squares sense by
    % using the backslash operator
    abc = [x y ones(length(x),1)] \ -(x.^2+y.^2);
    a = abc(1); b = abc(2); c = abc(3);

    % calculate the location of the center and the radius
    xc = -a/2;
    yc = -b/2;
    radius  =  sqrt((xc^2+yc^2)-c);

Alternative approaches will be appreciated...

Here is my solution, it is close to @Yoda's idea, but I changed some steps.

Mark all pixels such that there are at least 6 pixels in their 7x7 neighborhood
Remove all blobs, but the largest
Fill holes
Apply edge detection
Find circle using Hough transform

Here is the relevant Matlab code. I am using Hough transform for circles .m file in my code.

function FindCircle()
    close all;
    im = imread('C:\circle.png');
    im = im(:,:,2);

    ims = conv2(double(im), ones(7,7),'same');
    imbw = ims>6;
    figure;imshow(imbw);title('All pixels that there are at least 6 white pixels in their hood');

    props = regionprops(imbw,'Area','PixelIdxList','MajorAxisLength','MinorAxisLength');
    [~,indexOfMax] = max([props.Area]);
    approximateRadius =  props(indexOfMax).MajorAxisLength/2;

    largestBlobIndexes  = props(indexOfMax).PixelIdxList;
    bw = false(size(im));
    bw(largestBlobIndexes) = 1;
    bw = imfill(bw,'holes');
    figure;imshow(bw);title('Leaving only largest blob and filling holes');
    figure;imshow(edge(bw));title('Edge detection');

    radiuses = round ( (approximateRadius-5):0.5:(approximateRadius+5) );
    h = circle_hough(edge(bw), radiuses,'same');
    [~,maxIndex] = max(h(:));
    [i,j,k] = ind2sub(size(h), maxIndex);
    radius = radiuses(k);
    center.x = j;
    center.y = i;

    figure;imshow(edge(bw));imellipse(gca,[center.x-radius  center.y-radius 2*radius 2*radius]);
    title('Final solution (Shown on edge image)');

    figure;imshow(im);imellipse(gca,[center.x-radius  center.y-radius 2*radius 2*radius]);
    title('Final solution (Shown on initial image)');

end


Deconvolution of two FIRs Basic questions: What's the "correct" way to deconvolve two causal FIRs in the frequency domain (i.e. using the FFT), neither of which may be minimum phase but both may be considered to have stable exact inverses? By "correct", I mean, how does one avoid time-aliasing errors? If there is no "correct" way, what is the best way?

Some context: Given two causal FIRs, it is easy to avoid time-aliasing when convolving them in the frequency domain by first appropriately padding the IRs with trailing zeros (i.e. performing linear convolution via the FFT as discussed in standard DSP texts). However, it's not clear to me how to do something similar for deconvolution when the FIRs may not be minimum phase.

Deconvolving a minimum phase FIR out of a causal FIR (minimum phase or not) seems straightforward since we can compute the exact inverse of the minimum phase FIR (which we know will be causal and stable), and then linearly convolve this inverse with the other FIR.

However, since the exact inverse of a causal FIR that's not minimum phase is non-causal (assuming the inverse is BIBO stable - we force the inverse to have a finite length), performing linear convolution using the FFT after computing the inverse does not make sense since padding a non-causal inverse with trailing zeros doesn't make sense.

Related example with somewhat baffling results: In the figures below, I linearly convolve $x1$ with $x2i~(i = 1,2,3)$ to produce $yi$. Then, I deconvolve $x1$ out of each $yi$ to try and reproduce each $x2i$. To do this, I pad $x1$ with trailing zeros until the length matches that of $yi$, compute the exact inverse of the padded $x1$, and then circularly convolve the inverse with $yi$. The relevant difference between the plots across columns are the signals $x2i$ which have a different number of trailing zeros ($i$ denotes the number of trailing zeros). Can anyone provide some explanation as to why the results seem reasonable in columns 1 and 3 but not in column 2? Also, this might seem to be a reasonable way to deconvolve two signals and yet the results obtained don't always appear satisfactory (as this example shows). Any thoughts?

Note: In the plot titles, I have used $*$ to denote linear convolution and $/$ to denote deconvolution (all done in the frequency domain as described above).

Thanks!

I'm not sure why the problem occurs for you with the x22 example, but I believe the problem is because the deconvolution is not being performed over the right length.

Deconvolution of an FIR (inverting an FIR) requires, in general, convolution using an infinite duration impulse response.  So just using the default lengths in performing the deconvolution using an FFT might not generate the correct result.

I've implemented the same problem in R, but the last (additional) row of plots shows the results obtained if I do the deconvolution over a longer signal than the default. This generates the right results.

As you can see, my implementation generates bizarre results for the x21 example rather than the x22 example in yours.

R Code Below

x1 <- c(0,0,0,0,0,0,0,1,1,1,1,1,1,1,0,0)
x21 <- c(0,1,1,1,0)
x22 <- c(0,1,1,1,0,0)
x23 <- c(0,1,1,1,0,0,0)

my_conv <- function(x,y)
{
  out_length <- length(x)+length(y)-1
  x_pad <- c(x,rep(0,out_length-length(x)))
  y_pad <- c(y,rep(0,out_length-length(y)))

  z <- Re(fft(fft(x_pad) * fft(y_pad),inverse=TRUE))

  return(z/out_length)
}

my_deconv <- function(x,y, out_length=0)
{
  if (out_length == 0)
  {
    out_length <- length(x)+length(y)-1
  }
  x_pad <- c(x,rep(0,out_length-length(x)))
  y_pad <- c(y,rep(0,out_length-length(y)))

  z <- Re(fft(fft(x_pad) / fft(y_pad),inverse=TRUE))

  return(z/out_length)
}

par(mfrow=c(5,3))
plot(x1)
plot(x1)
plot(x1)
plot(x21)
plot(x22)
plot(x23)
y1 <- my_conv(x1,x21)
plot(y1)
y2 <- my_conv(x1,x22)
plot(y2)
y3 <- my_conv(x1,x23)
plot(y3)
plot(my_deconv(y1,x1))
plot(my_deconv(y2,x1))
plot(my_deconv(y3,x1))
plot(my_deconv(y1,x1,50))
plot(my_deconv(y2,x1,50))
plot(my_deconv(y3,x1,50))


Minimum sample frequency of IMU accelerometer and gyroscope I was wondering how I would justify a sampling rate of 125 Hz for accelerometer and gyroscope data from wearable sensors. This is a rate used in a lot of biomechanics literature, but I can't seem to find a justification.

I have been looking up the Nyquist rate, but I'm not sure how to find the bandwidth of the IMU signal.

Any help in this justification would be great

I'm not sure I can state anything better than Marcus's or Maxtron's comments (it's that way because it's always been that way / because the hardware is limited to that rate), but another way to think about it is:

The Nyquist rate is 62.5 (=125/2) Hz.
Can your acceleration or rotational angles (the things measured by the accelerometer and gyroscope) change at a frequency faster than that?

If your instrumented system cannot, then this sampling rate is OK. If your system can move faster, then the sampling will miss something.


Detecting objects in real world video feed? I have an image which consists on an object in a plain black background. This is the reference image of the object. Now, My task is to take a webcam and obtain the video feed and I have to detect if the object is there in the video feed.

For instance consider I have the image of a pen with a plain black background. Now when I take move my webcam around the office space I should be able to find the object.

What is the best approach for this?

There is no straight-forward answer to this. But I can tell you the things you need to consider to choose a method.
Before getting into the methods to detect the object in the reference image, let me tell you about invariance.

If the object in the video is of the same scale as reference image? ( I am guessing NO), then you have to consider detecting at various scale.

If the object in the video is of the same orientation as reference image? (I am guessing NO), then you have to choose a feature that is rotation invariant)

Now about various method that you can use. Remember this is non exhaustive and there are many other methods.

If you have only one reference image, then it is difficult to say without looking at the reference image and the kind of video you are referring to. But a simple answer would be template matching. Other method could be to compute feature descriptors like SIFT/SURF or even color of reference image and then try to search for it in the video. Or you can also use image processing techniques like thresholding, canny edge detection etc.

If you have 100s of images of same reference image, you can train a SVM classifier. The choice of feature depends on the object in the reference image.

If you have lots of different images(like in 1000s) of the reference image class, then try deep learning approach. If you have few few images, it is still possible with this dlib library for deep learning using MMOD(check this), but you need check the accuracy.

So it all depends on the problem in hand specifically.


2D Convolution as a Doubly Block Circulant Matrix Operating on a Vector I was reading Fundamental Image Processing, Chapter 5 (Image Transforms), I encountered the following problem:
Given the arrays $x_1(m,n)$ and $x_2(m,n)$ as follows:

Write their convolution $x_3(m,n) = x_2(m,n)\ast x_1(m,n)$ as a doubly block circulant matrix operating on a vector of size 16 and calculate the result.
 Can anyone please explain what the meaning of problem is and how it should be solved?

The point is that circular convolution of two 1-D discrete signals can be expressed as the product of a circulant matrix and the vector representation of the other signal.

The circulant matrix is a toeplitz matrix which is constructed by different circular shifts of a vector in different rows. For example, consider two signls $h[n]$ and $g[n]$, each of length $4$. If we assume $\mathbf{h}=\begin{bmatrix}
h_0 & h_1 & h_2 & h_3
\end{bmatrix}$, then the circulant matrix denoted by $\mathrm{circ}(h)$ is $$\mathrm{circ}(\mathbf{h})=\mathbf{H}=\begin{bmatrix}
h_0 & h_3 & h_2 & h_1\\
h_1 & h_0 & h_3 & h_2\\
h_2 & h_1 & h_0 & h_3\\
h_3 & h_2 & h_1 & h_0
\end{bmatrix}.$$
So we can calculate $h[n]*g[n]$ by evaluating the product $\mathbf{H}\mathbf{g}$, where $\mathbf{g}=\begin{bmatrix}
g_0 & g_1 & g_2 & g_3
\end{bmatrix}^T$.

This can be extended to 2-D signals. I explain it through the given example.

The size of signals are $2\times2$ and $3\times3$. The size of convolution is $(2+3-1)\times(2+3-1)=4\times4$. So we need to pad zeros to adjust the size of 2-D signal $x_2(m,n)$:
$$\mathbf{A}=\begin{bmatrix}
 0 & 0 & 0 & 0\\
0 & -1 & 0 & 0\\
-1 & 4 & -1 & 0\\
0 & -1 & 0 & 0
\end{bmatrix}$$
Now construct four circulant matrices corresponding to the four rows of $\mathbf{A}$ denoted by $\mathbf{a}_i$ as follows:

$$\begin{align}
\mathbf{X}_0&=\mathrm{circ}(\mathbf{a}_0)=\mathrm{circ}(\begin{bmatrix}
0 & -1 & 0 & 0
\end{bmatrix})\\
\mathbf{X}_1&=\mathrm{circ}(\mathbf{a}_1)=\mathrm{circ}(\begin{bmatrix}
-1 & 4 & -1 & 0
\end{bmatrix})\\
\mathbf{X}_2&=\mathrm{circ}(\mathbf{a}_2)=\mathrm{circ}(\begin{bmatrix}
0 & -1 & 0 & 0
\end{bmatrix})\\
\mathbf{X}_3&=\mathrm{circ}(\mathbf{a}_3)=\mathrm{circ}(\begin{bmatrix}
0 & 0 & 0 & 0
\end{bmatrix})
\end{align}$$
and construct the $16\times16$ doubly circuland matrix $$\mathbf{X}=\begin{bmatrix}
\mathbf{X}_0 & \mathbf{X}_3 & \mathbf{X}_2 & \mathbf{X}_1\\
\mathbf{X}_1 & \mathbf{X}_0 & \mathbf{X}_3 & \mathbf{X}_2\\
\mathbf{X}_2 & \mathbf{X}_1 & \mathbf{X}_0 & \mathbf{X}_3\\
\mathbf{X}_3 & \mathbf{X}_2 & \mathbf{X}_1 & \mathbf{X}_0
\end{bmatrix}.$$
Now we need to construct the vector corresponding to $x_1(m,n)$. Pad zeros from left and bottom to adjust indices considering $m,n$, and pad zeros from top and right to make the size $4\times4$:
$$\mathbf{B}=\begin{bmatrix}
 0 & 0 & 0 & 0\\
0 & 3 & 4 & 0\\
0 & 1 & 2 & 0\\
0 & 0 & 0 & 0
\end{bmatrix},$$
then construct a $16\times 1$ vector $\mathbf{b}=\begin{bmatrix}
 \mathbf{b}_0 \ |& \mathbf{b}_1 \ |& \mathbf{b}_2 \ |&\mathbf{b}_3
\end{bmatrix}^T$ by concatenating the rows of matrix $\mathbf{B}$.

Finally, the 2-D circular convolution is $$\boxed{\mathbf{G}=\mathbf{X}\mathbf{b}}$$


Matlab plot of QPSK system doesn't agree perfectly with theoretical BER curves Does anyone know if there is a simple explanation on the fact that the theoretical bit-error rate (BER) curves of a Quadrature phase-shift keying (QPSK) system are approximately 1 dB shifted from the simulated curves?

The simple explanation is that there is an error in your simulation. Here's one that works in MATLAB:

% number of symbols in simulation
Nsyms = 1e6;
% energy per symbol
Es = 1;
% energy per bit (2 bits/symbol for QPSK)
Eb = Es / 2;
% Eb/No values to simulate at, in dB
EbNo_dB = linspace(0, 10, 11);

% Eb/No values in linear scale
EbNo_lin = 10.^(EbNo_dB / 10);
% keep track of bit errors for each Eb/No point
bit_err = zeros(size(EbNo_lin));
for i=1:length(EbNo_lin)
    % generate source symbols
    syms = (1 - 2 * (randn(Nsyms,1) > 0)) + j * (1 - 2 * (randn(Nsyms, 1) > 0));
    % add noise
    syms_noisy = sqrt(Es/2) * syms + sqrt(Eb/(2*EbNo_lin(i))) * (randn(size(syms)) + j * randn(size(syms)));
    % recover symbols from each component (real and imaginary)
    syms_rec_r = sign(real(syms_noisy));
    syms_rec_i = sign(imag(syms_noisy));
    % count bit errors
    bit_err(i) = sum((syms_rec_r ~= real(syms)) + (syms_rec_i ~= imag(syms)));
end
% convert to bit error rate
bit_err = bit_err / (2 * Nsyms);

% calculate theoretical bit error rate, functionally equivalent to:
% bit_err_theo = qfunc(sqrt(2*EbNo_lin));
bit_err_theo = 0.5*erfc(sqrt(2*EbNo_lin)/sqrt(2));
figure;
semilogy(EbNo_dB, bit_err, 'bx', EbNo_dB, bit_err_theo, 'r', 'MarkerSize', 10, 'LineWidth', 2);
xlabel('Eb/No (dB)');
ylabel('Bit error rate');
title('QPSK bit error rate');
legend('Simulation','Theory');
grid on;

Note that the theoretical expression for bit error rate for BPSK/QPSK modulation is:

$$
P_b = Q\left(\sqrt{\frac{2E_b}{N_0}}\right)
$$

keeping in mind that $E_b$ is the energy per information bit. The somewhat subtle distinction between $E_b$ and $E_s$, the energy per symbol, is something that often trips up people new to the subject. This difference also explains why QPSK and BPSK have the same bit error rate when expressed as a function of $\frac{E_b}{N_0}$; you don't get any bit-error performance benefit by moving to QPSK, although you can achieve a given bit rate with less occupied bandwidth.


Output a sine wave from samples-limits? Using a micro-processor of 48MHz, I am trying to output a sine wave from the DAC (digital to analog ), and I am not sure if my Nyquist and calculations are good.

I have to output sine in frequencies range of 1KHz to 3MHz
I have created an array of a sine samples that has 1000 points(for 1 period )
I am running a timer between samples, and the frequency calculation

float timerFrequency=sineFrequency*1000; //(timerT*1000=finalSineT)

But then I was thinking that for 3MHz sine wave, I will need a timer frequency of 3GHz with this approach(!). (3Mhz signal * 1000 sine samples)

Something here is very wrong with everything, because if I have a processor that is 48MHz, I can probably output a 3MHz , isn't it ?

Whats wrong with my approach, and what are my limits with these assumptions ?
Can I reduce the sine samples array to lets say 100 ? what would it do to the Nyqvist at 3MHz ?
How many samples should my sine array needs? (the min), and my timer frequency?

Eventhough your problem is quite a practical one, based on a particular chip to produce a sine wave that requres a solid understanding of its blocks, I would, nevertheless, like to provide you a general framework for producing a continuous time sine wave from a set of discrete time samples, from which you can deduce the details of implementation for your particular architecture.

Now for this particular problem of producing a sine wave, the fundamental factor which will determine your analog output capability from a DAC is the quality of the analog "Anti-Imaging" filter (also known as interpolation or reconstruction filter)

When that filter is an ideal one (or a close enough approximation) it can reject all the distorting image frequency bands that exist due to the DAC output being an ideal impulse based sampled representation of the signal to be reconstructed. (A sine wave in this case).

In such a case (of close to ideal) of operating conditions, you can generate a perfect sine wave at any desired frequency from just "2" samples of it, by merely adjusting the output DAC conversion frequency. The filter will interpolate all the missing wave in the ideal case.

In other words, your sine wave samples will be this $$x[n] = \cos(\pi n)$$ which has only two samples which are +1 and -1 per period, representing a critically sampled sinusoidal. Your DAC will output +A and -A after scaling for each consequitive conversion period (which is limited by either your chip speed or your DAC output settling time) and that perfect reconstruction filter will convert those apexes into a smooth sine wave.

By adjusting the DAC output conversion period you can control the produced sinusoidal wave frequency simply as: $$f_{sine} = \frac{1}{2 T_{conversion}}$$

Now things are not ideal (especially DAC output involves a zero order hold block that should be taken into account in the design of the reconstruction filter) and therefore sadly just two samples per period will not be enough. In such a case, to be able to use a less perfect and easier to design analog output reconstruction filter you will provide more samples of the sine wave per period. But how many samples? It depends on your filter quality. The better the filter the less the number of samples.

The frequency is however is still adjusted the same way: In general let: $$x[n] = \cos(\omega_k n)$$
be your one period of sampled signal (look up table samples) where $w_k$ (being  between $0$ and $\pi$) is the discrete time frequency of the generic set of samples. And let your DAC output period be $T_s$ then the produced output sine wave will have the frequency of:
$$f_{sine} = \frac {\omega_k}{2 \pi T_s}$$

Again by adjusting $T_s$ or $\omega_k$ you can control this frequency.


g(x)=x odd and even expansions I'm studying book about digital filter by Richard Hamming. And there is exercise to get odd and even expansion of g(x)=x where x is from 0 to $\pi$. I understood even expansion, but can't get into odd expansion:
$\pi$ - $2$*(${\rm sin} x$ +  $\frac{1}{2}$${\rm sin} 2x$ + $\frac{1}{3}$${\rm sin} 3x$...)

Can anybody explain odd expansion of this function g(x)=x?

Thanks in advance.

UPDATE:
Let me be more clear. The exersise is: show that g(x) = x has two expansions (0 < x < $\pi$)

$x$ =$\begin{cases}\pi - 2*(\sin(x)+\sin(2x)/2+\sin(3x)/3+...)
\\\pi/2 - 4/\pi*(\cos(x)+\cos(3x)/3^2+\cos(5x)/5^2+...)\end{cases}$

I merely understand cosine expansion, but I don't understand sine expansion...

We want to find a representation in terms of sines and cosines for the function
$g(x)$ which has value $x$ for $x \in [0,\pi]$. We don't really care if this
representation does not work well outside the interval $[0,\pi]$, but
inside the interval, the representation better work perfectly.  This seems
to be crying out for a Fourier series approach, and so that's what we will
use.  The standard approach begins with functions that are periodic with
period $2\pi$ and so let us extend the definition of $g(x)$ to cover
one period $[0,2\pi]$.

If we choose $g(x)$ to have value $2\pi-x$ for
$x \in [\pi, 2\pi]$, then, when repeated periodically,
this extended $g(x)$ is an even function of $x$ with average value $\pi/2$.
Furthermore, as is the case with all even functions, the sine terms all
have $0$ coefficients.  Thus, the Fourier series works out to
have only cosines in it, and the DC term is $\pi/2$.
If we choose $g(x)$ to have value $x$ for $x \in [\pi, 2\pi]$, then
when repeated periodically,
this extended $g(x)$ is not an odd function of $x$, but nonetheless
it is true that all the cosine terms in the Fourier series
turn out to be $0$ except the DC term which has value $\pi$.


Is Periodic Convolution Same as Circular Convolution? I know exactly what is periodic convolution. But i don't know if circular convolution means to be the same thing!

You did not give specific details on how you define periodic or circular. Assuming standard definitions, the answer is yes.

Standard infinite-support data is treated with standard convolution. When the support  is finite, the "outer unknown samples" require assumptions.

Zero-valued outside samples is a common assertion. Another option, often used with Fourier or wavelet transforms, consists in arranging finite-support data on a ring or torus, the first sample following the last one.

It is often termed periodic, cyclic or circular, with a similar meaning.


Frequency response of numerical differentiation by polynomial interpolation / finite difference One can use polynomial interpolation (or finite difference) to do numerical differentiation. However, there seems to be a surprising lack of interest in obtaining frequency response of this numerical differentiation method. Is there any reference I can look up to, or can anyone provide what frequency response would be? I want to obtain frequency response of any filter order (any number of filter taps).

By higher order, I mean higher-order differences as in https://en.wikipedia.org/wiki/Finite_difference#Higher-order_differences or in case of polynomial interpolation, polynomial degree.

Polynomials are not integrable, hence their Fourier interpretation is complicated. However, locally, there are a few works. Under least-squares, polynomial interpolation's linearity works well with frequency analysis.

The class of Savitsky-Golay filters performs least-squares polynomial interpolation, and yields different derivative orders. As they seemed to be little known in signal processing, R. W. Schafer dedicated them What Is a Savitzky-Golay Filter? (public version) in the Lecture Notes section of IEEE Signal Processing Magazine in 2011. There, he especially covers their properties, including frequency responses. Analysis was also performed in Smoothing by Savitzky-Golay and Legendre filters by P. O. Persson and G. Strang in 2003. One can combine derivative accuracy and noise shaping or frequency response fitting.

More recently, there were derivations with IIR filters and fractional differentiation, for instance Maximally Flat IIR Smoothers With Repeated Poles and a Prescribed Delay, and even non-linear versions.

In image processing, there has been a lot of works as well, I can dig into my reference database if needed.


Difference between Analog IQ sampling and creating IQ digitally In the context of wideband receiver design, what are the advantages and disadvantages of Analog IQ sampling vs creating IQ digitally (with Hilbert transform)?

Why do we use different ADCs to create independent IQ channels when we can use Hilbert Transform to create the Q channel and use only one ADC?

You're right.  Whenever you can do things digitally you should because of all the advantages that you've listed.  But a long time ago it was expensive to do those low pass filters digitally.  If the LPF required is narrow (and requires a steep roll off), it would require many taps (multipliers) and they might require many bits so as to not lose precision to numerical accuracy issues. You also have to have the ADC twice as fast in the all digital case.
The disadvantages of doing the things in the analog realm are many though.  You would have to carefully tune LPFs to have the same group delay.  There could be many calibration issues.
Here's a paper describing a Hilbert Transform design


How/why are the $\mathcal Z$-transform and unit delays related? The $\mathcal Z$-transform uses the same notation as the unit delay $z^{-1}$, but in $\mathcal Z$-transform $z$ is a complex number.

What's the relation between the $\mathcal Z$-transform and the unit delay concept?
Also, is there a complex number interpretation of unit delay?

Mathematically we can easily show that the z-transform of a unit cycle delay is $z^{-1}$, Just like the s-transform of a time delay $\tau$ is $e^{-\tau s}$. I would like to add some additional insights that may lead to a more intuitive understanding:
Why z??
NOTE: To really understand how the z-transform and unit delays are related (intuitively), it really helps to understand the Laplace domain, what it is used for, and how it maps to the z-domain. For brevity I am including the first graphic covering a main point needed and the following paragraph giving a high level perspective, without going into further details of how we use either plane, but the more you know about that, the more the rest will make sense of what I will show after that.-- For further details on s and z along with a great graphic by Julius Smith that inspired my plots below, see  https://ccrma.stanford.edu/~jos/mdft/Comparing_Analog_Digital_Complex.html)
To help understand $z^{-1}$ it is useful to understand why we use the transform into the s and z plane. This is an entire topic on its own but I will at least make note of this very interesting observation: The process of multiplying (complex conjugate multiplication to be specific) and then integrating the product for two functions in the continuous time domain, or multiplying and then adding the product in the discrete time domain can be viewed as "correlation". (In actuality it is impressive how many things we know in DSP are a form of this description of correlation!).  Take a close look at the equations for the Fourier Series Expansion (FSE), the Fourier Transform (FT), Laplace Transform (LT) and Z-Tranform (ZT) equations and consider how we are doing correlations to each possible frequency for the FSE ($e^{jn\omega_nt}$), each point on the frequency axis for the FT ($e^{j\omega t}$), and each point on the complex surface for S ($e^{s t}$) and Z ($z^n$):
Fourier Series Coefficients
(scaled correlation of x(t) with $e^{jn\omega_0t}$, x(t) is only defined from 0 to T):
$$c_n = \frac{1}{T}\int_0^T x(t)e^{-jn\omega_0t}dt $$
Fourier Transform (correlation of x(t) with $e^{j\omega t}$) :
$$X(\omega) = \int_{t=0}^\infty x(t)e^{-j\omega t}dt$$
Laplace Transform (correlation of x(t) with $e^{st}$):
$$X(s) = \int_{t=0}^\infty x(t)e^{-st}dt$$
Z-Transform (correlation of x[n] with $z^n$):
$$X(z) = \sum_{n=0}^\infty x[n]z^{-n}$$
(x(t) and x[n] assumed to be causal, hence the integrations and summations start at 0).
(Note that for the Laplace and Z Transform it is not strictly a correlation the way I defined it since s and z have real components, so $e^{st}$ is not simply the complex conjugate of $e^{-st}$; both the real and imaginary exponent signs change.)
In each case, we have mapped our time domain function by describing it in terms of other variables, with the primary benefit of revealing its extrema (the poles and zeros) that can uniquely and compactly describe a (linear time invariant) system, opening up a whole world of elegant problem solving (notably solving integro-differential equations with simple algebra, and describing/ manipulating systems based on their pole and zero locations).
So we convert our discrete time domain sequence into the continuous z-domain by correlating it with every possible value of $z^n$ where z is any complex number, thereby decomposing our waveform into an infinite number of z's each with its own complex weight (magnitude and phase). This is the surface of our function on the z domain, which we usually just view the pole and zero locations (as that is sufficient to determine the magnitude and phase of every other point on that surface!).
The unit circle in the z-plane is the frequency axis:
Below demonstrates one way, out of many, to map from the s-plane to the z-plane. For purposes here, the only point I want to make is the reminder that the vertical axis in s which represents frequency $\omega$ (the Fourier transform is the function H(s) when you set $s=j\omega$),  maps to the unit circle $|z|=1$ on the z-plane.
(This mapping is specifically the matched-z transform, other common approaches are the bilinear transform, which I use most often, and the method of impulse invariance.  In all these mappings, what is consistent is the unit circle in z is the frequency axis (imaginary axis in s) and the left half plane in s maps to the inside of the unit circle in z. The answer on this link was very good for more details into the mapping techniques and why: Are there alternatives to the bilinear transform?)
Additionally this answer by Matt L describes additional mapping techniques where the $j\omega$ axis in the s-plane does NOT map to the unit circle in z-plane: z-Transform Methods: Definition vs. Integration Rule

Time Delay of One Sample = $z^{-1}$:
The math to show that the z-transform of a unit cylce delay is $z^{-1}$ is extremely simple (which is why we love the z-transform; we could certainly compute everything in discrete time systems using Laplace if we wanted the punishment).
$X(z) = Z\left\{x(n)\right\} =\sum_{n=0}^\infty x(n)z^{-n} $
$Z\left\{x(n-1)\right\} = X(z)z^{-1}$
So knowing this I provide the following graphics below to make the main point: as we increase powers of z, we are moving forward in time. The first graphic is of a shift register, where we often see each register denoted as $z^{-1}$, indicating that that element provides a 1 sample delay. Therefore the contents of the shift register at any given moment contain the current sample in the left most position, and all the previous samples as we move forward to the right (so we are moving backward in time reviewing the contents of our waveform in the shift register at that moment).  The second graphic shows the same waveform, but as we would be used to displaying it, on an axis that is moving forward in time, and that axis has increasing powers of z.

What is z???
Next, let's talk about z and what it is. For purposes here,the important point to know is that z is any complex number.

Now, knowing that z can be any complex number, and increasing powers of z represents samples as we move forward in time (and with the reminder when you multiply two complex numbers, you add the phases and multiply the magnitudes), let's look at $z^n$ as $n=0,1,2,3...$ for various values of z :

Notice in these plots how when z as any complex number is limited to the unit circle, the waveform that results for $z{^n}$ is a constant amplitude complex frequency, just as we introduced in the first plot how the unit circle represents the frequency axis.  When z=1, the frequency is DC, and when z is magnitude 1 with a very small angle, the frequency is small, and as the angle of z increases, the frequency increases. A delay in time of a complex number with magnitude 1 results in a rotation (the angles add), so if we continuously add delay with time, we will get a continuous rotation (frequency is change in phase vs change in time). This hopefully will give more insight in tying together an interpretive meaning for z, and $z^{-1}$ and how it relates to a unit sample delay.


Apply wavelet transform to analyse EEG signal I would like to apply The Morlet wavelet transform to analyse my EEG signals. I have many short signals each is only 1 min long. and they all recorded in 30Hz.
I have Two questions:

In the Morlet wavelet,  What is the best scale (alpha) to use in my case ?
About the edge effect: How can I know / compute what portion of my data will be corrupted due to the wavelet "edge effect" ?

I'm not sure whether you're talking about Discrete Wavelet Transform (DWT) or Continuous Wavelet Transform (CWT). Both can be used on discrete signals similarly to DFT and DTFT, I'm not sure if anyone calls it the Discrete Time Wavelet Transform instead. In any case, the dyadic wavelet transform is the non-redundant one (number of samples in $=$ number of coefficients out), and that's what we usually assume by talking about DWT.

It looks like you need to do a bit more reading on the wavelet transform in general, because of your first question. Wavelet transform doesn't work like the Fourier Transform; rather, it's akin to the Short Time Fourier Transform, which is a function of two variables: frequency and time. Similarly, the wavelet transform is a function of scale and time. This means that you don't pick a scale and stick to it. You go through many scales (in dyadic case they increment in powers of two), and calculate the transform for each one.

To address your second question, edge effects are usually not that easy to deal with and there's a plethora of papers on the topic. If you simply want to know what portion of the transformed signal will by affected by them, and it looks like that's what you're asking for, this paper has a good discussion. One thing to keep in mind is that the DWT, like the Discrete Fourier Transform, has circular wrapping symmetry, so dealing with edges is undesirable if you're hoping to get a perfect inverse transform as well. You can look at this paper for a discussion of that, and also a way to eliminate edge artifacts.


FM Modulation Index and Frequency trippler Suppose I have an FM signal with modulation index mf which is passed through a frequency tripler, then what will be the modulation index of output waveform ?.

Thanks in advance....

Hint: The modulation index is the peak frequency deviation divided by the maximum frequency of the baseband signal. Take a look at the formula for the FM signal and see what happens if the frequency is tripled. How does the peak frequency deviation change? How does the maximum frequency of the baseband signal change? I'm sure you can take it from here.

If not, please edit your question to let us know where you get stuck and ask a specific question (not just "please solve my problem").


Morphing between two FFTs I'm a musician and I'm making a script that takes a wave file and snaps each of its frequencies from the fourier transforms to the nearest musical harmonic. Thanks to help from another question I posted here, that part works, but what I need to do now is to experiment with some interesting ways to interpolate between the before and after version. I don't want to just do a crossfade because I want it to sound like one sound continuously morphing rather than one fading out and in. One idea I have is to snap one frequency at a time, maybe either in order based on the volume of that frequency by the end, from the lowest frequencies to the highest, highest to lowest, or random. Another idea was to interpolate them like image interpolation where if I have a frequency at 100hz with a volume of 255 and another at 120hz with a volume of 215, halfway between the interpolation I should have a frequency at 110hz with a volume of 235. If anyone knows of a function library that lets me do these then please let me know, otherwise I'll be making them and will need some help on getting started. How would I accomplish either of these?

Another thought I had was to reduce the number of harmonics over time so that in the beginning there are enough harmonics for the sound to be indistinguishable from the original but by the end there are only a couple. To get what I mean by the number of harmonics see the part of the script below where I calculate valid_ratios and valid_frequencies.

This is the sample I'm testing with (rename it to missile.wav):
https://my.mixtape.moe/iltlos.wav

Here is the script so far:

import numpy as np
import scipy.io.wavfile

samples, data = scipy.io.wavfile.read("missile.wav")

# fast fourier transform makes an array of the frequencies of sine waves that comprise the sound
data_fft = np.fft.rfft(data)

# the higher MAX_HARMONIC is, the more it sounds like the original,
# the lower it is, the more it sounds like an instrument
MAX_HARMONIC = 2

# generate list of ratios that can be used for tuning (not octave reduced)
valid_ratios = []
for i in range(1, MAX_HARMONIC + 1):
    for j in range(1, MAX_HARMONIC + 1):
        if i % 2 != 0 and j % 2 != 0:
            valid_ratios.append(i/float(j))
            valid_ratios.append(j/float(i))

# remove dupes
valid_ratios = list(set(valid_ratios))

# find all the frequencies with the valid ratios
valid_frequencies = []
multiple = 2
while(multiple < samples):
    multiple *= 2

    for ratio in valid_ratios:
        frequency = ratio * multiple

        if frequency < samples:
            valid_frequencies.append(frequency)

# remove dupes and sort and turn into a numpy array
valid_frequencies = np.sort(np.array(list(set(valid_frequencies))))

# bin the data_fft into the nearest valid frequency
valid_frequencies = valid_frequencies.astype(np.int64)
boundaries = np.concatenate([[0], np.round(np.sqrt(0.25 + valid_frequencies[:-1] * valid_frequencies[1:])).astype(np.int64)])
select = np.abs(data_fft) > 1
filtered_data_fft = np.zeros_like(data_fft)
filtered_data_fft[valid_frequencies] = np.add.reduceat(np.where(select, data_fft, 0), boundaries)

# do the inverse fourier transform to get a sound wave back
recovered_signal = np.fft.irfft(filtered_data_fft)

# amplify to ceiling
recovered_signal /= np.amax(recovered_signal)

scipy.io.wavfile.write("missile_output.wav", samples, recovered_signal)

maybe you can do it from bottom to top or top to bottom using something like the following.

a notation about notation:

$x[n]$ is the input sample for the $n$-th sample.  this originally comes from an A/D converter, but may have come from a sound file. $x_1[n]$ and $x_2[n]$ are two separate inputs.

$y[n]$ is the output sample at the $n$-th time.  it will eventually go to a D/A converter.

given subscripts $_1$ or $_2$ for inputs  $x_1[n]$ and $x_2[n]$, $a_k[n]$ is the amplitude of the $k$-th harmonic or the $k$-th sinusoidal partial (if the sound is not very periodic and has non-harmonic overtones or partials.  $\theta_k[n]$ is the instantaneous angle of the sinusoid that is the $k$-th partial, it's dimensionless but we can think of it as in units of "radians".  it is incremented at a rate of $\omega_k[n]$ (radians per sample), which is the instantaneous frequency of the $k$-th partial.

likewise for the output $y[n]$,  $b_k[n]$ is the amplitude of the $k$-th harmonic or the $k$-th sinusoidal partial, $\phi_k[n]$ is the instantaneous angle of the sinusoid that is the $k$-th partial, and $\nu_k[n]$ is the instantaneous frequency of the $k$-th partial of the output $y[n]$.

i am being lazy and not specifying the limits of $k$.  how many sinuosoidal harmonics or partials that you keep tabs on is a performance issue that is up to you.  too few and your sounds will be incomplete.  too many and you might have computational issues.  you don't want to treat every little spike or blip in the FFT spectrum as a bona fide "partial".  deal with the biggest peaks first.

for sound 1:

$$ x_1[n] = \sum_k a_{1,k}[n] \cos(\theta_{1,k}[n]) $$
where
$$ \theta_{1,k}[n] = \theta_{1,k}[n-1] + \omega_{1,k}[n] $$

and sound 2:

$$ x_2[n] = \sum_k a_{2,k}[n] \cos(\theta_{2,k}[n]) $$
where
$$ \theta_{2,k}[n] = \theta_{2,k}[n-1] + \omega_{2,k}[n] $$

and for the output:

$$ y[n] = \sum_k b_k[n] \cos(\phi_k[n]) $$
where
$$ \phi_k[n] = \phi_k[n-1] + \nu_k[n] $$

and let's say you have the harmonics ordered with $k$ the same for both inputs and output.  so you want $a_{1,k}[n]$ to morph to $a_{2,k}[n]$ and $\omega_{1,k}[n]$ to morph to $\omega_{2,k}[n]$ for all of the $k$.

so the amplitude envelope and frequency contours would have crossfade functions might look like:

$$b_k[n] \ = \ (1-g_k[n])\,a_{1,k}[n] \ + \ g_k[n] \,a_{2,k}[n]$$
or perhap you wanna do it to the logarithms of the amplitude, i dunno:
$$\log(b_k[n]) \ = \ (1-g_k[n])\,\log(a_{1,k}[n]) \ + \ g_k[n] \,\log(a_{2,k}[n])$$

similarly for frequency:

$$\nu_k[n] \ = \ (1-h_k[n])\,\omega_{1,k}[n] \ + \ h_k[n] \,\omega_{2,k}[n]$$

or log of frequency:

$$\log(\nu_k[n]) \ = \ (1-h_k[n])\,\log(\omega_{1,k}[n]) \ + \ h_k[n] \,\log(\omega_{2,k}[n])$$

so your question is, assumeably, how to define $g_k[n]$ and $h_k[n]$ so that the morph sounds like what you want it to.  is that correct?

now we'll introduce you to another simple (but conceptually salient) function that we audio guys call "clipping" or "limiting" but the computer graphics people call clamping:

$$ \operatorname{clamp}(x) \triangleq \begin{cases}
0 \qquad & x < 0   \\
x \qquad & 0 \le x < 1  \\
1 \qquad & 1 \le x  \\
\end{cases} $$

okay, so a possible way to morph this that might sound cool is to define the crossfading functions like this:

$$ g_k[n] =  \operatorname{clamp}\Big( \alpha(n-m) - \beta(k-1) \Big)$$

where $\alpha>0$ is your change rate and $\beta\ge0$ is the delay for the next harmonic to kick in..  $m$ is the sample index when the morph begins.  component number 1 will begin to change at sample $n=m$, but component number 2 will not change until sample $n=m+\beta$ and component 3 will not begin to change until time $n=m+2\beta$.  component $k$ will not begin to change until time $n=m+(k-1)\beta$.  do the same thing for $h_k[n]$.  this is the bottom frequencies changing before the higher frequencies.  you can make the top frequencies change first before the bottom frequencies by changing the sign on $\beta$.

you can also apply the Smoothstep function to the result of the $\operatorname{clamp}(\cdot)$ function.


Extremely computationally intensive image/video compression schemes which offer better compression than industry standard JPEG/H.264 I am trying to search for image/video compression schemes which offer superior compression performance (holding image quality as constant) when compared to widely adopted industry standards (i.e. JPEG/H.264), but for reasons of them being either too computationally or resource intensive, implementing them on CPU-based workstations is infeasible and therefore unadopted in common standards.

I could only think of fractal compression as a good example. Does the community know of any more of such examples?

There are various aspects of compression - whether you are referring to image or video. (May be audio is totally different so i am not referring this here).

If you really look at the history of compression standard, when they were first formed - MPEG1 real time encoders were rare. CPU at that time itself was not enough to make it real time; however, it was not the processing complexity (alone) that stopped more complexity algorithms to be used under codecs. (Probably the primary factor that was responsible for inclusion/exclusion of certain algorithms were a. patents, b. feasibility of implementation).

However, i would like to re-phrase the question a bit more simpler:

  What other research in compression is available which is still to come
  to market (or didn't come at all but might be worth for some
  applications)?

Here are some points to look at:

1. Better transform?

DCT by far has been universal transform across all codecs. However, there are many research papers which suggest that Wavelets could be better. Indeed the only JPEG2000 included wavelets - but that standard itself wasn't big hit. Probably it was better than DCT in some way; specially there are claims that Wavelets outperformce DCT [See: Ref 1]

I think wavelet is computationally expensive over DCT (partly because 8x8 DCT implementations are highly optimized). But i guess i didn't really flew that much because from the market perspective because DCT by the time was far more wide spread and the difference in improvement was not so substantial for people to switch.

2. Multi hypothesis motion compensation
While MPEG 2 was happy enough with $B$ pictures, extensive research from Bern Girod [See Ref 2] proved that - the concept can be extended way beyond this. Indeed H.264 actually has made great generalization on multiple target based motion compensaion - basically one can store past 16 pictures (compare to 2 in MPEG2) and indeed the prediction can be extremely high.

While, this is already part of the standard in H.264 - this is still not fully exploited. There is a need to make such extensively large search for motion estimation practically viable. Most H.264 encoder are still nascent in exploiting full potential of this.

3. Scalable Image/Video compression
Scalable Image and video compression is again a holy grail of how multiple description encoding can be done for image and video. Basically the same Image or video can be encoded at multiple Resolution, quantization (quality), or frame rate. These produces different layers of quality vs. bit rate right within the same stream - and the streaming server can actually adapt the bitrate and quality dynamically to best fit the video in the quality.
See [Ref 3], [Ref 4], [Ref 5]

Of this, only real success is progressive JPEG which is practically useful for websites and indeed useful. Again while the promise of dynamically adapting bit stream is great but it works only if there are practical systems where estimating such bit rate is possible and type of adaptation can serve some real purpose. However, in video domain, in spite the support of this in the standards (MPEG2, MPEG4 and H.264) - there is hardly any usage of this.

4. Object Video
This is perhaps ultimate holy grail of Compression. Object video [Ref 6] as initiated by MPEG standard body itself however, this was also backed by extensive research.

The basic idea is to segment the video in terms of distinct objects and later on each such object can be treated as separate plane. each object and be applied with different bit rate, quantization, frame rate and much efficient prediction can happen. Also, one can manipulate such objects like the way you can do stuff in photoshop/gimp.

However, as it turns out,  the segmentation of real-world image and video with perceptual relevant objects - is a hard problem! So even if i do have encoder/decoders - i need robus and real time segmentation system to make things possible. This is by far a really open problem as to how do you apply video in the form of such layered video to make it work.

5. More dimensions
Apart from this - there are codecs for stereo video (a.k.a multiview) and 3D video (almost similar concept) which additionally try to reduce the redundancy in the additional dimension. While this is still research going on here, H.264 does have a profile for this now. Same is true for the Cinema encoding.

There are many research projects still in the area of this - however, video compression is quite a commodity market by now. And research which doesn't substantially generate new application or if there are feasibility issues in implementation wont survive practically in spite it's ingenuity.

Hope this helps.

References:

A Comparative Study of DCT- and Wavelet-Based Image Coding by Zixiang Xiong, et. al. IEEE TRANSACTIONS ON CIRCUITS AND SYSTEMS FOR VIDEO TECHNOLOGY, VOL. 9, NO. 5, AUGUST 1999 pp. 692
Efficiency Analysis of Multihypothesis Motion-Compensated Prediction for Video Coding
Bernd Girod IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 9, NO. 2, FEBRUARY 2000 173
Spatial-domain resolution-scalable video coding by Barry G. Haskell, Hsueh-Ming Hang. Proc. SPIE 2094, Visual Communications and Image Processing November 718 (1993);
The MPEG-4 Fine-Grained Scalable Video Coding Method for Multimedia Streaming Over IP
by Hayder M. Radha, et. al. IEEE TRANSACTIONS ON MULTIMEDIA, VOL. 3, NO. 1, MARCH 2001 53
Advances in Scalable Video Coding by Ohm, J.-R.; Proceedings of the IEEE Jan. 2005 Volume: 93 Issue:1 pp. 42 - 56
Object-Based Texture Coding of Moving Video in MPEG-4, by Andr ́ Kaup, IEEE TRANSACTIONS ON CIRCUITS AND SYSTEMS FOR VIDEO TECHNOLOGY, VOL. 9, NO. 1, FEBRUARY 1999 pp. 5


Calculate the Output of Linear Time Invariant System Given it Impulse Response A filter is defined as $ h \left[ n \right] = \delta \left[ n \right] - \delta \left[ n - 1 \right] $.

Given a signal $ h \left[ n \right] $ defined as:

$$ x \left [ n \right ] = \begin{cases}
1 & \text{ if } x \geq 0 \\
0 & \text{ if } x < 0
\end{cases} $$

Let $ y \left[ n \right] = \left( x \ast h \right) \left[ n \right] $.
What is the value of $ y \left[ -1 \right], \, y \left[ 0 \right], \, y \left[ 1 \right], \, y \left[ 2 \right] $?

The Discrete Delta Function, $ \delta \left[ n \right] $ is the identity operator of Linear Time Invariant Systems.

Moreover, since it LTI System we can computer for each element of the filter by itself.
So the first element of the filter, $ \delta \left[ n \right] $, just outputs the signal itself.
The other element $ \delta \left[ n - 1 \right] $ just shifts the input signal.
Since the input is 1 for any $ n \geq 0 $ we subtract 1 from 1 unless it is $ n = 0 $ then we subtract zero from 1.

Hence the solution is 0, 1, 0, 0.


$\rm Q$ number format singed or unsigned? I am studying $\rm Q$ number formats and in a few places it's used to represent $2$'s complement fixed point numbers , and other places it's used to represent unsigned fixed point numbers. Given a $\rm Q$ number e.g $\rm Q1.4$:

How do I know if it's a signed or unsigned number ?
Should it be given that it's used in a signed or unsigned context?

From the name of the format alone ("Qm.n") you generally can't tell if the numbers are signed or unsigned. There is a convention that would use a prefix "U" to indicate that the numbers are unsigned, but that's not used universally. Often when the numbers are signed, the sign bit would not be included in the number m, so by comparing the register size with the sum m+n you could figure out if there's an extra sign bit or not. But that's also not a 100% indication because sometimes the sign bit is included in the number m, so m+n equals the register size, but the numbers are signed.

So in general you need extra information to be sure if the numbers are signed or not. In some cases, you can get a hint from comparing the register size with the sum m+n.


How do you plot an SPL vs frequency curve for a given set of pressure data? I have a set of pressure data at a sampling rate of 262144 Hz for 3 seconds making it 786432 samples of pressure data. I want to plot a curve of the SPL (sound pressure level) vs frequency. I first converted pressure data to SPL using $20\log_{10}(P/P_{ref})$ and performed a 4096 point FFT for the 192 segments of data and averaged them.

However, I couldn't get the required SPL values in decibels. Was the approach correct? If not, could you please point out the mistake in my approach?

The steps will be, in this order:

apply window function (e.g. Hann or Hamming)
calculate FFT
calculate magnitude of FFT (sqrt(re*re+im*im))
convert magnitude to dB (20*log10(magnitude))
apply combined correction/calibration factors (window function compensation, FFT factor of N scaling if needed, factor of 2 for symmetric FFT, Pref calibration, etc) - this can just be a single dB correction value that you add to the dB magnitude values

If your Pref calibration data varies with frequency (e.g. microphone with non-flat frequency response) then it's a little more complicated than the above, but the same basic principles apply.

Note also that if you want weighted dB SPL (e.g. A weighting) then you'll also need to apply the relevant frequency dependent weighting corrections.


Any idea what's the OCR tool behind `myfonts.com/WhatTheFont` font recognizer? I went to myfonts.com/WhatTheFont to get the font from a text to train tesseract.
I realized that the font recognizer did actually better OCR than when I used tesseract on the same image.
So any idea how could I find out what OCR tool or library they're using?
They do have a forum but it's not a support forum, it's more like a font marketplace.

I believe this question is off-topic, but this page suggests that the algorithm is described in:


Banding in the derivative of a lock-in amplifier signal I'm analysing several measurements taken with a SR830 lockin amplifier.

These measurements look similar to this one.
Since I'm interrested in the derivative of the signal I took the numerical difference. What I found was this.
One can clearly see two bands in the signal. I'm no expert, but I'm guessing it might be related to digitalisation in the lock-in amplifier.

So this leads me to two questions:

What could be the origin of these two bands and how to avoid them in the first place?
How should I correct the data already taken?

It might not be so straight forward to simply take the derivative of the output from a lock-in amplifier. If you think about how a lock-in works, the final stage of a lock-in amplifier is the low-pass filter. So if you then simply take the sample-to-sample derivative, you will simply be undoing the work of the low-pass filter and getting noise. Instead, you need to decimate your time domain signal, and then take the derivative.


how to eliminate interference in the same center frequency I'm trying to receive an ASK signal in the 433 MHz center frequency with a bandwidth of 100KHz. I'm using a CC1101 transceiver connected to a Raspberry Pi GPIO ports. However, there is an intermittent signal in the 433MHz center frequency coming from the sensors of my alarm system, so the RX doesn't get all the samples of my wanted signal and I'm not sure how I could filter the interference because both signals are "unknown" (they are not generated by me). I wonder how I could do this, at least in theory.

First of all, the CC1101 supports a number of other frequencies, so as a trivial solution, you can switch to using those.

The other thing that you can do at the RF module level is to "shift" your data out of the problematic region. Now, I am assuming that the other module is on 433MHz and it is not a similar module with similar operation. Let's assume that the interfering module is different than the one you are using.
What you could do, is shift to 4-FSK with the highest DEVIATION_M (refer to the datasheet) but use it as 2-FSK. What you would be trying to do is use the two symbols that cause the largest deviation to the centre carrier frequency to avoid your data being contaminated by the interference. Usually, these two symbols would be 00 and 11 (notice these pairs). This is inefficient of course, but it might work as a quick workaround for you.

But what might be happening is that the other RF module is also behaving in the same way. It might also be doing Frequency Shift Keying, in which case, none of the frequencies your RF module might try to hop to is "safe" from noise.

In this case, you have to work with high levels of noise.

The number one way to guard against this kind of interference is with the use of coding.

This can take many forms and have various complexities. One of the simplest "encoding" schemes is a Repetition Code coupled with simple Majority Decoding scheme. You might notice above that your effective B-FSK used the symbols 00 and 11 which, effectively, is a bit like using a repetition code.

One of the issues you might face in this case is that you do not have access to the internals of the RF module. In other words, the only thing you would be doing here is repeat bits $N$ times and at the reception end, receive a frame and decide what bit was sent by assigning the frame to the majority of bits received (either $0$ or $1$).

Obviously the question now is, where does the frame start and where does it end?

To deal with this, you are now going to have to embed your data stream on some form of protocol that specifies that your actual message starts with something like 0000000000111111111100001111001101 followed by your encoded bits.

The frame start sequence has equal chances of being contaminated with noise as the rest of the data, but the idea is to be looking in to your incoming data stream for a sequence that looks as close as possible to the frame start sequence.

And thus, progressively, we are going closer and closer to convolutional codes where you would be inserting redundant 0s and 1s in a predictable way in your data stream to be able to decode the most likely sequence to have been sent. That is a bit more involved at the decoding side but convolutional codes are more powerful than plain repetition codes.

Hope this helps.


Do symmetric discrete signals have zero phase? I generated a Hanning window having an even symmetry:

where for even-sampled case we either take "left" and "right" to include or exclude the center sample. I was surprised to find its DFT has non-zero imaginary coefficients, i.e. non-zero phase, in both cases. I've also tried an odd symmetry case and ended up with non-zero real coefficients.
Shouldn't even- & odd-symmetric discrete signals have zero-phase DFT?

Assuming you are referring to the standard definitions of odd and even symmetry, i.e.
\begin{align}
\text{even: }&& x[-n] &= x[n] \\
\text{odd: }&&x[-n] &= -x[n]
\end{align}
then it still holds that the DFT of an even sequence is real (i.e. the phase is 0 or $\pi$) and the DFT of an odd sequence is imaginary (i.e. the phase is $\pi/2$ or $-\pi/2$).
This is independent of whether the the length of the sequence is even or odd as long as the sequence meets the definition of symmetry. This implies that for odd sequences $x[0]$ must be $0$ and for odd sequences of even length $x[N/2]$ must be zero too, where $N$ is the DFT length.
Both sequences that you show in your original question are not "symmetric" in the mathematical sense. They are only symmetric when properly time shifted and time shifting is clearly a non zero-phase operation.
The following Matlab script illustrates that this is true within the numerical accuracy of double precision math:
%% script to confirm zero phase DFT of even/odd symmetric sequences

N = 128; % FFT length

%% Even length sequences

% create a random sequence of length N/2+1
x = randn(N/2+1, 1);

% create even version
xEven = [x; x(end-1:-1:2)];

%  perform DFT and check zero phase
fx = fft(xEven);
zeroPhaseError = 10 * log10(sum(imag(fx).^2) / sum(abs(fx).^2));
fprintf('Even symmetry, Even length: zero phase error = %6.2fdB\n', zeroPhaseError);

%% Odd length sequences

% create odd symmetry version. This requires the samples at n=0 and n=N/2 to be
% zero
xOdd = [0; x(2:end-1); 0 ; -x(end-1:-1:2)];
%  check zero phase
fx = fft(xOdd);
zeroPhaseError = 10 * log10(sum(real(fx).^2) / sum(abs(fx).^2));
fprintf('Odd Symmetry, Even Length: imaginary phase error = %6.2fdB\n', zeroPhaseError);

%% now try an odd length sequence
% create a random sequence of length N/2
x = randn(N/2, 1);

% create even symmetry version
xEven = [x; x(end:-1:2)];

%  check zero phase
fx = fft(xEven);
zeroPhaseError = 10 * log10(sum(imag(fx).^2) / sum(abs(fx).^2));
fprintf('Even symetry, odd length: zero phase error = %6.2fdB\n', zeroPhaseError);

% create odd symmetry version. This requires the samples at n=0 to be zero.
% There is no sample at n=N/2 since that's not an integer.
xOdd = [0; x(2:end); -x(end:-1:2)];
%  check zero phase
fx = fft(xOdd);
zeroPhaseError = 10 * log10(sum(real(fx).^2) / sum(abs(fx).^2));
fprintf('Odd symmetry, odd length: imaginary phase error = %6.2fdB\n', zeroPhaseError);


Identification of the system with zeros I need to identify two systems. Those are their Unit Step responses:

I don't know, what order of the system that is. It could be second order system with zero in 0. Lecturer told me it was third order system, but then told me to solve it as second order system, so I'm a bit confused. More over we don't have materials for those systems, which would show us how to idnetify them. Lecturer told me to search it on the internet.

So I wonder whether anyone could give me some hint or better, if anyone could explain how to solve this or what kind of system those systems are.

I would welcome some good materials for this task.

Hint: It's a system with a damped oscillation of a single frequency sine wave, which also starts with an arbitrary phase. That should give you amplitude curve (damping), oscillation component (sinusoidal wave) and phase (starting point of the wave). If you can figure out the general form this system takes mathematically, it's a matter of plugging in the right values for those three.

Another hint:
Don't overthink the system order. Think about what a first order system could look like. Do you need to go further? What could a second order system look like? Do you need to go further, etc. I guarantee you will stop pretty early.

Perhaps a more important hint:
You're given a step response, not an impulse response, but there's a way to go from one to the other (what is it). Once you find the expression to describe this curve, going to impulse response should be very easy.

One last hint:
Depending on how your subject is taught, this may help or this may hinder your attempts, so take this one with a grain of salt. If it confuses you, ignore it. My Linear Systems professor said this and made us repeat this out loud in class almost every lecture: Solutions to linear systems are always exponentials!


Is it possible to get the velocity vector for tangential movement using radar? When an object appears towards the radar sensor or goes away from the radar sensor, there is a change in the phase shift of the frequency and from that, the velocity can be inferred.

So, I have two questions.

Is it also possible to find out the velocity "vector" ( ie which direction the object is moving to and not just the magnitude of the velocity of the object )
If the object is neither moving towards the sensor nor away from the sensor ( moving tagentially ) , will there be still a phase shift on the radiated pulse?

The answer to your first question is yes if you can resolve the change in bearing or azimuth as well as range in a sequence of measurements, and some variation of a Kalman Filter is used.

The answer to your second question is a lot more interesting because Special Relativity predicts a tangential Doppler shift.  One common scenario used in presenting Special Relativity to students, involves two rocket ships flying at different near light speed velocities comparing their clocks.

This may interest you.

H. W. Thim, "Absence of the relativistic transverse Doppler shift at microwave frequencies," in IEEE Transactions on Instrumentation and Measurement, vol. 52, no. 5, pp. 1660-1664, Oct. 2003.
doi: 10.1109/TIM.2003.817916
Abstract: An experiment is described showing that a 33-GHz microwave signal received by rotating antennas is not exhibiting the frequency shift ("transverse Doppler effect") predicted by the relativistic Doppler formula. The sensitivity of the apparatus used has been tested to be sufficient for detecting frequency shifts as small as 10-3 Hz which corresponds to the value of (v/c)2 = 5.10-14 used in the transverse Doppler shift experiment reported here. From the observed absence of the transverse Doppler shift it is concluded that either the time dilation predicted by the standard theory of special relativity does not exist in reality or, if it does, is a phenomenon which does not depend on relative velocities but may be a function of absolute velocities in the fundamental frame of the isotropic microwave background radiation.
keywords: {Doppler shift;microwave antennas;microwave measurement;33 GHz;absolute velocities;isotropic microwave background radiation;microwave frequencies;relativistic transverse Doppler shift;rotating antennas;time dilation;Doppler effect;Doppler shift;Masers;Microwave antennas;Microwave frequencies;Microwave measurements;Microwave theory and techniques;Receiving antennas;Testing;Transmitting antennas},
URL: http://ieeexplore.ieee.org/stamp/stamp.jsp?tp=&arnumber=1240186&isnumber=27804


DFT calculations returning different results My coworker ran a DFT query using a software in C downloaded to his PC. His results was the following chart:

I am using a JS tool from corbanbrook/dsp.js

My results were very different.

any ideas why, or where to look?

Make sure you use the same scaling and X and Y axis so you compare apples to apples. The upper graph is in power the lower seems to be amplitude
Start with a known signal with a known answer. A sine wave or square wave for example and compare the outcome with the expected result.


Understanding 2D alignment using least squares In Computer Vision Algorithms and Applications book by Richard Szeliski, the author explains 2D alignment using least squares as follows:

  Given a set of matched feature points ${(x_i, x_i')}$ and a planar parametric transformation of the form
  $$x' = f(x;p),$$
  how can we produce the best estimate of the motion parameters $p$? The usual way to do this
  is to use least squares, i.e., to minimize the sum of squared residuals
  $$E_{LS} = \sum{\lVert{r_i}\rVert^2} = \sum{\lVert{f(x_i,p)-x'}\rVert^2},$$
  where $r_i$ is the residual between the measured location $x_i'$ and its corresponding current predicted location $f(x_i; p)$.

  Many of the motion models presented in Section 2.1.2 and Table 2.1, i.e., translation,
  similarity, and affine, have a linear relationship between the amount of motion $\Delta{x} = x'-x$ and the unknown parameters $p$,
  $$\Delta{x}=x'-x=J(x)p$$
  where $J=\partial{f}/\partial{p}$ is the Jacobian of the transformation $f$ with respect to the motion parameters $p$ (see Table 6.1). In this case, a simple linear regression (linear least squares problem)
  can be formulated as
  $$E_{LSS} = \sum{\lVert{J(x_i)p - \Delta{x_i}}\rVert^2}$$

Question:
In the last equation, what is the point of trying to minimize the difference between $J(x_i)p$ and $\Delta{x_i}$? What is a higher-level interpretation?

The set of feature points $\boldsymbol{x}_i$ and $\boldsymbol{x}^{\prime}_i$ is given, and you want to optimize the parameters $\boldsymbol{p}$ of a given planar transformation $\boldsymbol{f}(\boldsymbol{x};\boldsymbol{p})$ such that it predicts the actual transformation from $\boldsymbol{x}_i$ to $\boldsymbol{x}^{\prime}_i$ with minimum error. In this way, the parametric transformation optimally describes (and predicts) the actual transformation.

One way to do this is to minimize the error between the actual difference $\Delta\boldsymbol{x}_i=\boldsymbol{x}^{\prime}_i-\boldsymbol{x}_i$ and the difference predicted by the chosen transformation. The predicted difference is formulated in terms of the Jacobian of the transformation $\boldsymbol{f}(\boldsymbol{x};\boldsymbol{p})$ in order to get a linear least squares problem (which is easy to solve): ${\Delta\widehat{\boldsymbol{x}}_i}=\boldsymbol{J}(\boldsymbol{x}_i)\boldsymbol{p}$.

So the error to be minimized is given by

$$\sum_i\left[{\Delta\widehat{\boldsymbol{x}}_i}-\Delta\boldsymbol{x}_i \right]^2=\sum_i\left[\boldsymbol{J}(\boldsymbol{x}_i)\boldsymbol{p}-\Delta\boldsymbol{x}_i \right]^2\tag{1}$$

Chosing the parameters $\boldsymbol{p}$ that minimize $(1)$ minimizes the average squared error between the actual feature point difference $\Delta\boldsymbol{x}_i$ and the one predicted by the transformation $\boldsymbol{f}(\boldsymbol{x};\boldsymbol{p})$.


is there a way to apply a median filter diagonally? For example, some studies show how to smooth a spectrogram with a 2D median filter vertically and horizontally (along time or frequency). My question is how would I apply 2D median filter along the diagonals of the spectrogram instead? This would be useful for a certain study I am doing, so I need to understand how that would be coded in MATLAB.

The median filter can be applied on many pixel-set shapes, and indeed this was used to emulate of fasten full 2D median filters by a combination of leaner medians.

In Matlab, your request just amounts to gather diagonal elements in a vector. The simplest could be:

extract each sliding square patch of the image
extract the diagonal as a matrix with diag.m
compute the median

If you want the antidiagonal, flip the square patch with flipud, fliprl.
If you want other diagonals than the $\pm 45°$, you have to use discrete lines, for instance with the Bresenham algorithm.

If you want to remain high level, I'd suggest you to use gblk: a pedestrian data block extractor at MatlabCentral to parse an image with overlapping blocks.


Designing a discrete-time inverse system I am trying to design a discrete-time inverse system to eliminate an undesired echo in a data transmission problem.

The echo caused by the transmission channel is represented as attenuation by a factor of 0.9 and a delay corresponding to one time unit of the input sequence.

The received signal $y[n]$ can be expressed as follow;

$$y[n] = x[n] + 0.9x[n-1]$$

So basically, I'm trying to determine the unit impulse response of a causal inverse system to discover $x[n]$ from $y[n]$.

Conceptually, I'm thinking I need to connect a system $h_2[n]$ in series with the system $h_1[n]$ such that $h_1[n]*h_2[n] = \delta[n]$.

Can anyone offer some advice on how best to go about this ? I've been searching online but can't find a concrete example to use as a starting point.

Correct answer is in the comment already but for completion: Your transfer function is

$$H(z) = b_0 + b_1 \cdot z^{-1}, \mbox{ with } b_0 = 1, b_1 = .9$$
which makes the inverse
$$H^{-1}(z) = \frac{1}{b_0 + b_1 \cdot z^{-1}} $$
which corresponds to the difference equation
$$y(n) + 0.9 \cdot y(n-1) = x(n) $$ or
$$y(n) = x(n) -0.9 \cdot y(n-1)  $$


Best way to detect steps and slopes in a signal what is the best way of detecting larger steps and slopes in a signal?
(as shown in the pictures, 1. to be detected signal, 2. good signal)
I thought about some kind of sliding window analysis.
Thanks for help.

This problem is heavily related to the notion of piece-wise constant or piece-wise linear approximations. Either to fit the experimental curve at best, or to find the optimal parameters may requires some more details, like: do you know the number of parts? Any idea about the noise? Some literature to start with:

Best piecewise constant approximation of a function of single variable
Models and Algorithms for Optimal Piecewise-Linear Function Approximation
Fitting Piecewise Linear Continuous Functions


Viterbi Decoder - Decoding Continuously Encoded Stream Partially I am in trouble with Viterbi decoder for a while. I have some questions.

I need to decode messages from a continuous stream which is encoded without resetting the state metrics. Each message starts with a preamble, but since these preambles are also encoded as part of the stream, there is actually no message concept before decoding.

Encoding without resetting the metrics means that the data is encoded as one big message, doesn't it?
The decoder I am using works frame by frame, without storing and re-using the state metrics. So, I am adding a prefix and a postfix from the previous and next messages.

(last $n$ bits of the previous encoded message + message to be decoded + first $m$ bits of the next encoded message)

I hopefully expect that the prefix works as a training data, and the state metrics after the prefix are ready to decode the actual message. In order to decode the last bits of the message correctly, the postfix is used. I suppose using the data from the next message as the postfix is better than adding zeros or random bits, because adding zeros or random numbers may cause a wrong path affecting not only the decoding of the postfix but also of the actual message.

Is this procedure reasonable?
Is it common and does it have a name such as sliding windows?

How many bits should I use as the prefix and postfix? (I cannot set the traceback length and see if it is pre-set somewhere. The code rate is $1/2$ and constraint length is $7$.)

Thanks in advance.

Since the decoder that you are working with is a black box, yes, adding prefixes and postfixes to the stream to "window-ize" it is appropriate.  The pre and postfixes need to be the actual data from the previous and post windows, and they need to be at least as long as the traceback depth.  If you do this and get rid of the result bits that correspond to the prefixes and postfixes, then you should do fine.

In a normal situation you would not need to add any postfixes or prefixes.  There are only two things you need to do to handle continuous streams of Viterbi-encoded data.

1) occasionally check the state metrics and lower them when all of them get "high".  This works because the only thing that matters is the difference between the state metrics, not their absolute value.  An easy way to do this is to check if all of them are equal to or above half the maximum state metric value.  If they are, subtract half the maximum value from all of them.

2) You have to do tracebacks while you are decoding.  This is do-able in both FPGA's and microprocessors (it's actually quite easy in microprocessors).  Every $N$ transitions (pick an $N$ that works for your architecture) you do a traceback and output $N$ bits.  You have to go back through $M$ bits ($M$ being equal to the traceback depth) before you get to the $N$ bits that you will output.


Force an unstable IIR filter to be stable by forcing a nonlinearity in the digital block diagram I am modeling an analog filter with digital software and have reduced the model to a 4th order FIR filter in discrete space with transfer function

$$ H(z) = \frac{b_0 + b_1 z^{-1} + b_2 z^{-2} + b_3 z^{-3} + b_4 z^{-4}}{1 + a_1 z^{-1} + a_2 z^{-2} + a_3 z^{-3} + a_4 z^{-4}} $$

It happens to be unstable, but I have noticed that most analog "musical" LP filters distort non-linearly when in this realm of high resonance rather than blow up to infinity (of course).
Is it possible to simulate this by adding a nonlinear element to the IIR block diagram?

Clipping with $y[n] = \operatorname{clip}(x[n] + \ldots)$ creates a periodic [-1, -1, 1, 1, ...] in the output signal, so is there another way I can prevent numerical blowup?

You won't be able to avoid the instability with a linear model. With your current discrete model, the filter will blow up at a certain resonance setting. Using a better filter topology or higher precision will allow you to go a bit farther, but it will blow up.

Preventing that would require a saturating nonlinearity in the feedback path. That path is hidden in your implementation. To access it, you need to go back to the 4-section model without feedback and add the nonlinearity when you add the feedback path. But then you won't be able to simplify the filter to a fourth order linear filter.

Your options are then to delay the feedback path with its nonlinearity by one sample and feed it back into the 4-section model input. That will create significant deviations from the intended filter behaviour for high frequencies and high resonances. You can compensate for that somewhat, but it is tedious to get right. See Antti Huovilainen's paper for an approach like that.

The other option is to solve the implicit feedback equation using an iterative method for every single sample. That is doable, but requires significantly more computational effort and some experience with numerical methods to guarantee convergence and stability.

If you go that latter route, you might as well go with a proper nonlinear filter model, because that will need very similar methods to solve.


Segmentation and Labeling of White Blood Cells I am trying to segment and label white blood cells from this image:

Using K-Means I managed to segment most of what I wanted:

But as you can see, there are some cells connected to each other which I would like to label as separate objects.

So my questions is:
Can you point me to some known algorithms that perform the task of automatically separating them?
This image comes from a data set with annotations of the cell centers, in case machine learning is involved...

Also, but much less important, it would be nice to hear suggestions about improving the K-Means segmentation (Discarding the 'Stains').

Edit:
Thanks!
After playing around with the regional minima function
(inextendedmin), I managed to segment all of the cells in the center, except for the bottom left in that cluster at first, but some cleaning up of the segmented image did the trick in the second watershed segmentation:

Edit2:
For anyone who still follows (or even cares), I managed to improve segmentation using
K-means for spatial clustering on the original bw image, where the number of cells I obtained via watershed was defined as k (number of clusters). The result provided a bit finer segmentation:

The cells image is taken from ALL-IDB

I can't recommend anything in the realm of your clustering algorithm, but this looks to be pretty simple. The initial segmentation worked nicely; a simple watershed should separate those cells.

For example, in MATLAB:

function W = dowatershed(I)
    %% I is a logical mask (i.e., im2bw)
    d = -bwdist(~I);
    l = watershed(d);
    W = I;
    W(l == 0) = 0;
end


Multirate Signal processing for matching DAC of the SDR I have been breaking my head for a long time. Please help with this..any input would be helpful!

I am trying to implement a OFDM system using SDR (USRPN210 for Transmitter and USRP2 for receiver) using Simulink and my laptop (core i5).

Baseband specifications:
OFDM signal sample rate = 25.6MS/s
No. of samples in one frame = 10880
duration of one ofdm frame = 425e-6s
signal bandwidth = 25.6MHz
Used bandwidth is 22.4MHz

USRP specs:
DAC Master clock rate (fixed)= 100MHz
Interpolation factor = ??

As far as I have understood, I have to design a multirate system which converts 25.6MS/s to a sample rate which is a factor of the master clock rate i.e. 100MHz. And give the factor as the Interpolation factor.

For example, convert 25.6MS/s to 25MS/s and give interpolation factor as 4, so it matches with the 100MHz(25*4). However this causes Underrun. So I need to increase the interpolation factor and decrease the input baseband sample rate.

My problem is how to design a downsampler which will reduce the 25.6MS/s to a lower sample rate without loss of data.?

25.6 MHz is larger than the largest possible sampling rate of these two USRP models (25MHz). Therefore, this is impossible to implement directly; in any way, you'll need to generate your 25.6 MS/s signal first, and then resample to 25  MS/s. That resampler will be very CPU-consuming, probably more than your whole OFDM signal processing.

The underruns happen because your computer is, including that resampler, simply not fast enough to compute the samples at that rate. You'll need a faster computer, a more efficient resampler, or a different OFDM configuration.

Now, I don't know which resampler you're using, but if it's one of these included in GNU Radio, there's only two options:

the PFB arbitrary resampler with a ratio of 0.9765625
the rational resampler with an interpolation of 125, and a decimation of 128.

Both resamplers are really CPU-hungry, but for your case, I'd actually expect the arbitrary resampler to work slightly more efficient than the rational resampler in their default configurations, where the filters are choosen automatically based on the rate ratio alone.

However, you have a signal with unused carriers at the Nyquist band edges - you don't need as steep a core filter in your classical multirate resampler as dictated by the resampling ratio. So, make a drawing of the aliases that a decimation filter would have to suppress, and design a filter that works there. It will be significantly shorter than what would be designed automatically!


Feature extraction for sound classification I'm trying to extract features from a sound file and classify the sound as belonging to a particular category (eg : dog bark, vehicle engine e.t.c). I'd like some clarity on the following things :

1) Is this doable at all? There are programs that can recognize speech, and differentiate between different types of dog bark. But is it possible to have a program that can receive a sound sample and just say what kind of a sound it is? (Assume there's a database containing a lot of sound samples to refer to). The input sound samples can be a bit noisy (microphone input).

2) I assume that the first step is audio feature extraction. This article suggests extracting MFCCs and feeding them to a machine learning algorithm. Is MFCC enough? Are there any other features that are generally used for sound classification?

Thank you for your time.

By long shot it is doable - to what extend? You will see. This task of environmental sound classification is not very well studied. Also choice of machine learning paradigm is crucial - statistical approach or maybe binary classifier? You can start with GMM's, ANN's and SVM's - I opt for GMM's and ANN's.
Yes, most of people are using MFCC's because they are well correlated with what people are actually hearing and also no one came up with anything better since. You might also want to add extra features such as MPEG-7 descriptors.  Proper feature optimisation must be performed because sometimes you don't need so many features, especially when they are do not separable. For more info please refer to my previous answers:

Feature extraction from spectrum
MFCC extraction
Detection of sounds


Discrete Time Signal Property Property:

The frequency of oscillation of discrete time sinusoids sequence increases as $\omega$ increases from $0$ to $\pi$. If $\omega$ is increased from  $\pi$ to $2\pi$ then frequency of oscillation decreases.

My question:

What is meant by "frequency of oscillation" ?

How can a frequency of a signal vary (increase) from $0$ to $\pi$ and (decrease) from $\pi$ to $2\pi$. Isn't it should be same from $0$ to $2\pi$ ?

I think it means the "apparent" frequency of oscillation. It's poorly worded.

What's happening is that the frequency is increased to the Nyquist rate in the first half and then above it, causing aliasing, in the second.

What looks like the frequency increases in the first half and decreases in the second because of aliasing.


SNR relation with modulation index and noise power What is the relation between the AM modulation index and the SNR? (I plotted this and it came as the shape of the exponential graph)

Also what is the relation between the channel noise power and the SNR?
I plotted this and  I got the shape as this graph:

(ignore the labels of the axis)

If $x(t)$ is your message signal, then you can write the AM signal as

$$s(t)=A[1+mx(t)]\cos\omega_0 t$$

where $A>0$ is a real-valued constant, $m$ is the modulation index, and $\omega_0$ is the carrier frequency (in radians). If you assume that the noise added by the channel is white with power spectral density $N_0$, then the SNR after demodulation is

$$SNR=\frac{A^2m^2\overline{x^2(t)}}{4N_0B_x}$$

where $B_x$ is the bandwidth of the message $x(t)$ and $\overline{x^2(t)}$ is its average power.


What are good ways to detect signal clipping in a recording? Given a recording I need to detect whether any clipping has occurred.

Can I safely conclude there was clipping if any (one) sample reaches the maximum sample value, or should I look for a series of subsequent samples at the maximum level?

The recording may be taken from 16 or 24-bit A/D converters, and are converted to floating point values ranging from $-1...1$. If this conversion takes the form of a division by $2^{15}-1$ or $2^{23}-1$, then presumably the negative peaks can be somewhat lower than -1, and samples with the value -1 are not clipped?

Obviously one can always create a signal specifically to defeat the clipping detection algorithm, but I'm looking at recordings of speech, music, sine waves or pink/white noise.

I was in the middle of typing an answer pretty much exactly like Yoda's.  He's is probably the most reliable but, I'll proposed a different solution so you have some options.

If you take a histogram of your signal, you will more than likely a bell or triangle like shape depending on the signal type. Clean signals will tend to follow this pattern. Many recording studios add a "loudness" effect that causes a little bump near the top, but it is still somewhat smooth looking. Here is an example from a real song from a major musician:

Here is the histogram of signal that Yoda gives in his answer:

And now the case of their being clipping:

This method can be fooled at times, but it is at least something to throw in your tool bag for situations that the FFT method doesn't seem to be working for you or is too many computations for your environment.