Improve coverage of BackgroundSubtractorMOG2 and findContours I'm trying to track players on a soccer field, using a fixed camera. The picture below shows the results of BackgroundSubtractorMOG2, the second picture below the results of findContours on the results of BackgroundSubtractorMOG2.

As you can see there are a lot off small contours detected. How can a I alter the results of BackgroundSubtractorMOG2 so one shape (a player) results in one contour I can track ?

This is the code I used:

import numpy as np
import cv2

cap = cv2.VideoCapture("images/Keeper.mov")
history = 30   # or whatever you want it to be
accelerate = 5

kernel = cv2.getStructuringElement(cv2.MORPH_ELLIPSE,(3,3))
fgbg = cv2.BackgroundSubtractorMOG2()

while(1):
    for i in (1, accelerate):
        ret, frame = cap.read()

    fgmask = fgbg.apply(frame, learningRate=1.0/history)

    fgmask = cv2.morphologyEx(fgmask, cv2.MORPH_OPEN, kernel)
    cv2.imshow('frame',fgmask)

    h, w = fgmask.shape[:2]

    contours0, hierarchy = cv2.findContours( fgmask.copy(), cv2.RETR_TREE,
    cv2.CHAIN_APPROX_SIMPLE)
    contours = [cv2.approxPolyDP(cnt, 3, True) for cnt in contours0]

    vis = np.zeros((h, w, 3), np.uint8)
    levels = 17
    cv2.drawContours( vis, contours, (-1, 3)[levels <= 0], (128,255,255),
            -1, cv2.CV_AA, hierarchy, abs(levels) )
    cv2.imshow('contours', vis)

    k = cv2.waitKey(30) & 0xff
    if k == 27:
        break

cap.release()
cv2.destroyAllWindows()

I'm not an expert in this field, but have you tried eroding and then dilating?
I'm writing a python application to track tennis players in the court, and what I did in order to get a single contour for the players is to erode and dilate the tracked shapes.

erosion = cv2.erode(img,None,iterations = 1);
dilation = cv2.dilate(erosion,None,iterations = 1);
cv2.findContours( dilation.copy(), cv2.RETR_TREE, cv2.CHAIN_APPROX_SIMPLE);

or you can directly call the function

opening = cv2.morphologyEx(img, cv2.MORPH_OPEN, kernel);
cv2.findContours( opening.copy(), cv2.RETR_TREE, cv2.CHAIN_APPROX_SIMPLE);

Source: OpenCv Morphological Transformations


Can compressed sensing be used instead of intepolation for missing values? Consider a signal that is sparse in frequency, but it measured in the time domain, for example (in matlab):

t=linspace(0,4,2^8);
x=sin(2*pi*t)+sin(2*pi*3*sqrt(2)*t+pi/4);

let's say that some fraction of the signal is missing, for example the 120-130 measurement:

x(120:130)= NaN;

using standard interpolation methods (nearest neighbors, spline, etc) will fail to reproduce the missing signal, because these methods do not make use of frequency information of the signal.

I thought that by using compressed sensing (L1-magic for example) and random sampling of the signal outside the missing part I'll be able to recover the signal and the missing part, but looking into it I saw that I still need the full signal to create a measurement vector via y=Ax, where A is the sensing matrix.

If x has NaN values then y will be all NaNs, and the entire process is compromised.

So how can I "compress" my sampling by not having the full signal even though it has the information needed? (assuming the sparsity above)...

I'm aware I may use other techniques but I'm curious how to use CS.

Yes, at least in the above case it is possible.
Though it might not be computationally as cheap as other methods such as least squares based curve fitting methods.
I do not think injecting NaN gonna help, instead  let's look at the ignored data as dimension reduction feature of CS.
Here, the measurement matrix is basically an identity matrix that some of its row (120:130) have been removed:
A = eye(256);
A(120:130) = [];

Just be aware, for the case of measurement matrix you posed in your question, it might violate coherence requirement.
If we apply the above matrix, we cut away samples of 120 to 130. Now assume the resultant as measurement vector $y$ and $A$ as measurement matrix. This is how CS recovery algorithms work, e.g. in case of Fourier sparse signals, we have an incomplete set of measurements $y$ ,the goal is to find sparsest set $a$ such that if we apply A matrix, yields something close enough to $y$.
N = 2^8;
t=linspace(0,16,N);
x=sin(2*pi*t)+sin(2*pi*3*sqrt(2)*t+pi/4);
x=x';
Phi = eye(256);
Phi(10:17,:)=[];
y = Phi * x;
Psi=dctmtx(N)';  % Transposed DCT matrix is like inverted DCT
recovered=SolveOMP(Phi*Psi,y); % any solver available!
CS_intepolate=idct(recovered);
plot(CS_intepolate,'*'); hold on; plot(x,'r')

(blue line is CS interpolated and red is original signal)
BTW, as far as I understood, you suggest to mix the remaining data. See it as system of equations, blending the current equations into new ones, does not add any new information or constraint to the system, so it is not helping.
Also I believe there is connection with CS and data completion. You might find it interesting also.


Does the duration of a signal affect its frequency component's amplitude? Also, does the sampling frequency affect the power of a signal? I have two questions that are bugging me:

Does the duration of an audio signal affect the amplitude of frequency components of that same signal? For example, I am recording the sound of a fan using a microphone. At first, I record only for 10 sec and convert the audio signal into frequency spectrum. Then, I record the same sound for 20 sec and then convert the audio signal into frequency spectrum. In both the cases, the sound of the fan is same, but does the duration of the signal affect the amplitude of frequency components in the spectrum plot?
For example, I have 2 audio signals. For the first one, I have that same fan sound recording for 10 sec and the sampling frequency is 5KHz, and for the second recording, I have that same audio signal but now the sampling frequency is changed to 15KHz. I used MATLAB to check the power for both the signals and the power for both the signals was same, however I want to know why. Formula that I used was  Power=rms(signal)^2. According to me the second signal should have more power because now there are more number of samples compared to the first recording and since those extra samples would also have a random amplitude, the average shouldn't be the same as for the first one. Am I thinking it right?
Can anyone provide their thoughts? Thank You!

Power is energy per unit time. If you increase the duration, you increase the energy, but due to the normalization with time the power would be the same.
The DFT as given by
$$X[k] = \sum_{n=0}^{N-1}x[n]e^{-j 2\pi n k/N}$$
will scale the frequency component by $N$ as given by the summation over $N$ samples. This can be normalized by multiplying the result by $1/N$.
The frequency components of the signal levels will be the same in a normalized DFT (normalized by dividing by the total number of samples) for signal components that occupy one bin (pure tones), but the noise floor as observed may be lower by the change in sampling rate: if the noise floor is limited by quantization noise, the total quantization noise (well approximated as white noise, meaning constant across all frequencies) will be spread over a wider frequency range, so increasing the sampling rate will cause the contribution of quantization noise on a per Hz basis (noise density) to be lower. Further the duration will effect the frequency resolution in each bin in the frequency domain; for an unwindowed DFT, the equivalent noise bandwidth per bin is $f_s/N$ where $f_s$ is the sampling rate and $N$ is the number of bins.  At a given sampling rate, increasing the number of bins will increase the total duration and thus reduce the noise bandwidth per bin; causing the overall DFT noise floor as observed to decrease. (Same noise power just less measured in each bin). Windowing if done in the time domain will reduce the signal level by the coherent gain of the window, but increase the equivalent noise bandwidth such that pure tones will be reduced more than noise that is spread over multiple bins. This detail is explained further in this post.


Why use transfer functions than differential equations? I have some simple questions regarding DSP and Control Systems, but found no simple answers. I just need simple and easy to understand answers, not complex answers with huge maths.

Why we use Transfer Functions, when we can get a system's output by just solving it's differential equation?
What we finally get after solving a transfer function (In Laplace Transform)?

Why we use Transfer Functions, when we can get a system's output by just solving it's differential equation?

Because differential equations are unwieldy and hard to deal with, and you can't see the behaviour on different frequencies from these, whereas transfer functions just give you the behaviour of an LTI system given an excitation of given property. Much easier to work with!


Bandwith of a energy signal I'm studying for my finals , and I'm struggling with finding bandwidths.

I have this energy spectral density : Ψ(f)=rect(f) , and I need to find the bandwith, 3 dB band and the band with 90 % power.

About the 90 % case, do I need to find the total energy and then ∫W0Ψ(f)=0.9×TotalEnergy?

If I had a power signal and wanted to find the 3 dB bandwidth, would I need to find the average power and then solve ∫W0Sx(f)=0.5×TotalPower ?

Thanks in advance.

The concept of the bandwidth has several definitions.

For simple filters, one definition of bandwidth is the -3 dB frequency (for a lowpass filter) $f_c$ where the filter gain at $f_c$ is $H(f_c) = H(0)/\sqrt(2)$; i.e. $f_c$ is the half-power or cut-off frequency. This definition could be applied for a bandpass filters too but not for highpass filters. Note again that this definition of bandwidth requires that the filter frequency response be a continuous waveform, otherwise the the frequency will be undefined.

For signals, a 3dB bandwidth concept will not work and instead something like a 90 % energy (or power) definition is invoked. In this case the frequency spectrum of the signal is used to find the minimum frequency $f_b$ at which 90% of the total enegy of the (baseband, real) signal is contained.

In your case you have a signal with a rectangular frequency spectrum, which defies a 3 dB bandwidth, but enables the computation of a 90 % bandwidth definition.

However note that all these definitions make sense under typical waveform assumptions, but for your case neither -3dB nor 90% enegy definitions adequately define the signal's bandwidth.

You signal is bandlimited in the sense that the ractangular function is exactly zero out of its lobe and that frequency shall be taken as the bandwidth of this signal.


How to evaluate a background subtraction algorithm given ground truths? I'm testing a background subtraction algorithm with sample image sequences. I am wondering how one can evaluate the accuracy of a background subtraction algorithm results given ground truths?

The only idea I have in mind right now is to take the difference between a result image and its corresponding ground truth image and calculate difference/total_white_pixels_of(ground_truth).

Cheers!

EDIT:
Dear all, by ground truth image I mean binary image with foreground as white blobs and background as black pixels. This is an example of input image and its ground truth:

source: http://wordpress-jodoin.dmi.usherb.ca/dataset2014/

One way would be to create a simple signal-to-noise estimate, using the pixel intensity and the ground truth image as a mask. Here's an example of what I mean in Python:

import numpy as np

# Sample data
# Indices [0, 1, 2, 6] are background, [3, 4, 5] are object
original  = np.array([1, 2, 2, 5, 6, 8, 2])
truth     = np.array([0, 0, 0, 1, 1, 1, 0])
back_mask = 1 - truth
result    = np.array([0, 1, 0, 4, 5, 7, 1]) # Result of background subtraction

def non_masked_mean(input_array, mask):
    """mean of non-masked elements in the array"""
    return np.ma.masked_array(input_array, mask).mean()

def snr(input_array, back_mask):
    """ratio of non-masked to masked mean intensities"""
    with np.errstate(all='ignore'):
        snr = non_masked_mean(input_array, back_mask) / non_masked_mean(input_array, (1-back_mask))
    return 0 if np.isnan(snr) else snr

snr_before = snr(original, back_mask)  # (mean of objects / mean of background), before
snr_after  = snr(result, back_mask)    # (mean of objects / mean of background), after
snr_ratio  = snr_after / snr_before

if snr_after == np.inf:
    print "Perfect SNR in result! Either great or suspect..."
elif snr_after == 0:
    print "Image was flattened :("
else:
    print "SNR changed by a multiple of %.2f" % (snr_ratio)

The problem with this is that you get a perfect score by favouring a heavy-handed background subtraction. Even if only 1 pixel of the object survives this SNR will be infinite. Another approach would be to consider the background and the signal separately, and weight the results of the two. Then you can decide if it's more important to flatten the background, or preserve the original.

def non_masked_mean_ratio(input_array, result_array, mask):
    """ratio of two array means, with the same mask"""
    return non_masked_mean(result_array, mask) / non_masked_mean(input_array, mask)

signal_ratio     = non_masked_mean_ratio(original, result, back_mask)     # (mean of objects after) / (mean of objects before)
background_ratio = non_masked_mean_ratio(original, result, truth)         # (mean of background after) / (mean of background before)
weight           = 0.5   # Higher weight favours preservation of objects

score = (weight*signal_ratio) + ((1-weight) * (1-background_ratio))

This gives the score a nice expected range of about 0-1, where 1 is the ideal result. It's not perfect though, as you can boost it by increasing the intensity of the objects.

If you're expecting a binary result, you can use a similar strategy except using the number of non-zero pixels instead of the mean:

def nonzero_count_ratio(result_array, ideal_array):
    """ratio of actual number of non-zero pixels vs ideal non-zero pixels"""
    return float(np.count_nonzero(result_array*ideal_array)) / np.count_nonzero(ideal_array)

true_obj_ratio   = nonzero_count_ratio(result, truth)     # Number of correct object pixels vs number of potential correct
false_back_ratio = nonzero_count_ratio(result, back_mask) # Number of false background pixels vs number of potential false
weight           = 0.5

score = (weight*true_obj_ratio) + ((1-weight) * (1-false_back_ratio))


FFT Algorithm: Split Radix vs Radix-4 I wanted to know the difference between split radix and mixed radix algorithm.

I also wanted to know what is the difference between radix 2 ,radix 4 and radix 8 FFT

just did a Google image search and these diagrams have popped up.

these are radix-2 butterflies with expanded detail (showing where the "twiddle factor", Q, goes) for various DIF and DIT FFTs.

a radix-2 FFT (the simplest in concept) uses only radix-2 butterfly operations, each with 2 complex numbers going and and 2 complex numbers coming out.

these four radix-2 butterflies combined

make a radix-4 butterfly:

the image on the right just above would be a short-hand notation for a radix-4 butterfly.  if there were 8 signals going into the dot and 8 coming out, it would be a radix-8 butterfly.

let $p \in \{1, 2, ...\}$ and $P \triangleq 2^p$.  essentially a radix-$P$ FFT is an FFT that uses radix-$P$ butterflies and a radix-$P$ butterfly is, in and of itself, a complete $P$-point FFT with some phase adjustment on the internal twiddle factors.

a split-radix FFT has butterflies of different sizes.  here is the best example i can think of:

a 256-point FFT can be done entirely using radix-4 butterflies because 256 is a power of 4: $256=4^4$.  there would be 4 FFT passes and all of the passes would have radix-4 butterflies.  a 256-point FFT can also be done with a radix-2 FFT and it would be 8 passes, all with radix-2 butterflies, instead of 4 passes.

now suppose you wanted a 512-point FFT.  if you started with a radix-2 FFT, all you need to do is double the size of the buffer and add one more pass because $512=2^9$.  so there is no split radix if you wanted to do that.  turns out that the radix-4 FFT is a little more efficient than the radix-2 so you might want to replace 8 contiguous passes of the radix-2 FFT with 4 passes of a radix-4 FFT and one final (or one initial) pass using 2-point butterflies.  that is a split-radix FFT.


Is up-sampling prior to cross-correlation useless? Consider a simple case where two signals from two different sensors are cross-correlated, and the time-delay-of-arrival computed from the absissa of the peak of their cross-correlation function.

Now let us further assume that due to the dimensionality constraints of both antennas and the constraints on maximum possible sampling rate, the maximum attainable delay possible is $D$, corresponding to 10 samples.

The problem:

Because of those constraints, your computed delay may vary from any integer value between 0 and 10 samples, that is: $0 \le D \le 10$. This is problematic because what I really want is fractional-delay discrimination of the delay between the two signals impinging on my antennas, and changing the dimensions or the sampling rate are not an option.

Some thoughts:

Naturally, the first thing I think of for this case is upsampling the signals before performing a cross-correlation. However I think this is 'cheating' somehow, because I am not really adding any new information into the system.
I do not understand how upsampling is not 'cheating' in a sense. Yes, we are reconstructing our signal based on its currently observed frequency information, but how does this give one knowledge of where a signal truly started between, say, $D=7$ and $D=8$? Where was this information contained in the original signal that determined that the true fractional-delay start of the signal was actually at $D=7.751$?

The question(s):

Is this truly 'cheating'?

If not, then where is this new 'information' coming from?
If yes, then what other options are available for estimating fractional-delay times?

I am aware of upsampling the result of the cross-correlation, in an attempt to garner sub-sample answers to the delay, but is this too not also a form of 'cheating'? Why is it different from upsampling prior to the cross-correlation?

If indeed it is the case that the upsampling is not 'cheating', then why would we ever need to increase our sampling rate? (Isnt having a higher sampling rate always better in a sense than interpolating a low sampled signal?)

It would seem then that we could just sample at a very low rate and interpolate as much as we want. Would this then not make increasing the sample rate 'useless' in light of simply interpolating a signal to our heart's desire? I realize that interpolation takes computational time and simply starting with a higher sample rate would not, but is that then the only reason?

Thanks.

It's not cheating, and it's also not adding any new information.  What you are doing is the same thing that any upsampling LPF is doing- adding zeros and then reconstructing the waveform with the already known frequency information.  Thus, there is no new information, but there is still finer time resolution.

Upsampling the result is similar- no new information but finer time resolution.  You can do something very similar through quadratic interpolation.

All of these methods- upsampling and polynomial interpolation- get their information on where the fractional peak is from both the peak itself and its neighbors.  A quick pictorial example.

The blue line in the picture above is my simulated cross-correlation data (though it could be any result, not just a cross-correlation).  It is what I call a "balanced" peak because the neighbors are symmetric.  As you might expect, the resulting quadratic interpolation (red line) indicates that the true peak is at zero.

The image below, on the other hand, shows an unbalanced peak.  Please note that nothing has changed in the result except for the values of the two nearest neighbors.  This causes the interpolator, though, to shift its estimate of the fractional peak.

A nifty side benefit of these methods (polynomial interpolation and upsampling) is that it also gives you an estimate of the true peak value, though we are usually more interested in the location.

  If indeed it is the case that the upsampling is not 'cheating', then
  why would we ever need to increase our sampling rate?

To satisfy the Nyquist criterion.

  Isn't having a higher sampling rate always better in a sense than
  interpolating a low sampled signal?

No.  From a theoretical standpoint, as long as the Nyquist criterion is satisfied it doesn't matter what the sample rate is.  From a practical standpoint you generally go with as low a sample rate as you can get away with to reduce the storage requirements and computational load, which in turn reduces the resources that are needed and power consumption.


Bose Chaudhuri Hocquenghem (BCH) vs Reed Solomon (RS) Error correcting codes and sound transmission channel modeling I've implemented an acoustic based information transmission system using a bandwidth of 2 kHz and the band between 18.5 and 20.5 kHz. I'm trying to find the best error correcting code for it, in the sense that I could get the same error correction capability with the lowest redundancy added (I want to get the best speed possible). I have been able to transmit data between 2 smartphones at low rate (400 bps) and short distances (5-10 cm). At the moment I've design a package that I send between peers and I'm using 4 BCH(63,39,4) blocks there. I have also added some parity check bits (just for error detection) in case BCH fixed the incoming package in the wrong way (could happen if more that 4 errors occur in any of my 4 BCH blocks). I interleave the BCH coded bits before I modulate them trying to equally distribute the error burst among the 4 BCH blocks. Right now I'm capable of implementing both RS and BCH coding schemes, but I was wondering if there was some criteria I could follow in advance for choosing which one might perform better, would it be worth the time implementing RS. I read this document which says that for a Rayleigh Fading Channel being simulated in Matlab, BCH outperformed binary RS for a specific N,K,T selection. But they mention that this has a strong channel model dependency. My questions are:

What channel model fits data transmission over air best, I guess distortion generated by speakers should also be take under account, specially at the frequencies at which I'm working?
Is there a general study comparing BCH vs RS in different channels?
Is RS theoretically better than BCH for the same code rate n/k. Non binary RS may be considered, I'm just interested in same error correcting capability with greatest bit-rate ?

Independent answers to any of these questions are most welcome

Note: At this low bitrate CPU is not a bottleneck. Anyway the most time consuming task of my system is FFT based matched filter demodulation (specially if I use high order modulations), not BCH error decoding with such small code length.

What channel model fits data transmission over air best, I guess distortion generated by speakers should also be take under account?

Disclaimer: I'm really not an Audio person. There's a lot of audio DSP engineers on here that probably have a far better idea of what's happening on an acoustic channel. But:

Frequency selectiveness will be pretty important. It's unlikely your phone speakers can go below 40 Hz, or far above 12 kHz.
Audio typically might experience a bad case of multipath
prepare yourself for what would be a relatively solid amount of Doppler if your communication partners, or significant reflectors, or the medium, are moving

  Is there a general study comparing BCH vs RS in different channels?

probably, hundreds. But: it's really a thing that needs to be compared under specific channel models, and since there are freely available implementations of both, as soon as you have a realistic channel model for your specific use case down, simulation should give you more application-specific results than any study might ever give you.

I'd really go for channel sounding, in your case; send a known signal of suitable bandwidth, record the raw reception, try to extract a transfer function from that (assume LTI, and verify that assumption). Change the environment/channel, and do it again. Do that ten times (and write a paper), and come up with a stochastic channel model (e.g. "we have this general amplitude response for the ideal single channel, we have x~(PDF) multipath components, and their delays and weights follow the XYZ distribution", and write a paper).

Then write a channel simulator (at this point, I might, again, recommend GNU Radio, for which several channel sounders and channel simulators exist, but mainly for RF channels – but that's not inherent to its design, just to most of its userbase; obviously, write a paper), and test a few codes on that (and write a paper :))!


detecting "deformed" lines I have a project in which I have to recognize 2 lines in different "position", the lines are orthogonal but can be projected on different surfaces. I'm using opencv.

I'm not asking for  code, but only for hints about how can I solve this? I tried houghlines function but it works only for straight surfaces.
thanks in advance !

The cross shape that i see in this example seems to remain more or less orthogonal. Maybe using a houghlines transformation with a piecewise flexibility might provide more true positives, but this also increases the false positives. In this case a corner detector might also be a goodidea, thus providing a set of interest points to search from and eliminate. The problem again there would be the large number of false positives. Try this morphological corner detector (Robert Laganiere (1998): "Morphological Corner Detection") - its simple and structural. The problem in such cases is classifying your background on which the cross falls to help improve with the false positives!


python: how to compute the gray level histogram features as mentioned in the paper, and Hi I am extracting the grey level features of image mentioned in this paper (part 4.5 low level features).

  "We describe two low level features that are particu- larly important for photo quality assessment – contrast and brightness. Professional photos usually have higher contrast than snapshots. We measure the contrast of an image as fol- lows. First, we compute the gray level histogram Hr , Hg , and Hb for each of the red, green, and blue channels, re- spectively. Then, we compute the combined histogram H, where..."

I can't clearly understand the above part in the paper.
As my understanding about grey level features, take the above image as example.
Use the opencv convert the color image to grey image. as the following code:

image = cv2.imread('lenna.png')
gray = cv2.cvtColor(image, cv2.COLOR_BGR2GRAY)

and then using the following code to get the histogram with the following code:

  hist = cv2.calcHist([gray], [0], None, [256], [0,255])

I want to ask, and I don't know the meaning of gray level histogram for each channel.

How to we compute the gray level histogram Hr , Hg , and Hb for each of the red, green, and blue channels?
How to compute the the contrast quality, qct mentioned in the paper?

(1) how to we compute the gray level histogram Hr , Hg , and Hb for each of the red, green, and blue channels.

The gray level histogram for each of the channels is simply the gray level histogram of the Red channel, the Green channel and the Blue channel separately. Instead of converting the image from RGB to grayscale, try to extract its RGB components.

For example, in GNU Octave, a colour image is represented as a $M \times N \times D$ matrix, where $M$ is the height of the image, $N$ is the width of the image and $D$ is the "number of channels".

In a colour RGB image, the number of channels is 3.

This can be viewed as 3 images packed together in one "3D" matrix. To work on the Red channel, all that you have to do is use slicing to obtain just one of those matrices.

In Octave, this would be done via:

image = imread("someimage.png")
red_channel = image(:,:,1)

Similar thinking applies to other languages / computing platforms.

  2) how to compute the the contrast quality, qct mentioned in the paper.

According to the paper, the $q_{ct}$ is the width of the 98% "part" of the combined histogram.

So:

Load the image
Get the histogram of the Red Channel (call it $H_r$)
Get the histogram of the Green Channel (call it $H_g$)
Get the histogram of the Blue Channel (call it $H_b$)
Produce the combined histogram as $H = H_r + H_b + H_g$
Get the 98% of its mass:

Sort $H$ in reverse order of values (largest count comes first)
Get the cumulative sum of the histogram values until you have covered the 98% of the total sum.
Find the minimum and maximum values of the indices of the values that contribute to the sum
Calculate the $q_{ct}$ as the difference between the maximum and minimum of the indices.

Hope this helps.


How to implement Polyphase filter? I am trying to implement the following 4-Tap Polyphase window shown in figure below and used in : polyphase filters- radio astronomy: lecture #8 by Prof. Dale E. Gary on Matlab. I have tried splitting the window in time domain after multiplication with Sinc function but has not achieved the sidelobe level and mainlobe width shown in green graph- figure on right. I am not sure what I might be missing here.

I assume that the Polyphase filter shown in the paper is used as a decimator such that the output rate is an integer fraction the input rate (and this would make sense in a correlation operation since the output frequency required is less due to the proper filtering of the correlator. Making a polyphase filter implementation is quite easy; given the desired coefficients for a simple FIR filter, you distribute those same coefficients in "row to column" format into the separate polyphase FIR components as explained in the following example:

Assume a FIR filter with 8 taps for simplicity of explanation, with  coefficients as follows:

h[0], h[1], h[2], h[3], h[4], h[5], h[6], h[7], ...

and we use that to create a 4 element polyphase each with 2 taps:

Poly FIR1: h[0],   h[4]

Poly FIR2: h[1],   h[5]

Poly FIR3: h[2],   h[6]

Poly FIR4: h[3],   h[7]

From the above it is hopefully clearer in what I meant by "row to column" format.

The Polyphase Decimator would be structured as in the figure below:

A Polyphase Interpolator works the same way, and since I already have the graphics that explain clearly why we do the row to column mapping, I will include that below (and the beauty of the polyphase implementation itself intuitively explained!). This shows an input signal of successive pulses at a higher rate going through an interpolation operation. Normally, without the use of polyphase implementations, we can interpolate a signal by simply inserting zeros, and then following that with a low pass filter to get rid of the higher frequency aliases that occur due to the zero insert. The best filter design would pass our original spectrum with no distortion and completely eliminate the higher frequency components. It is the coefficients of this FIR filter that are used in the polyphase implementation since it mathematically performs the same result (that said I suspect the filter chosen in the paper based on windowing method of filter design is not the best choice).


Why are magnitudes normalised during synthesis (IDFT), not analysis (DFT)? In most examples and FFT code that I've seen, the output (frequency magnitudes) of the forward DFT operation is scaled by N -- i.e. instead of giving you the magnitude of each frequency bin, it gives you N times the magnitude.

Operationally, this is simply because the DFT is calculated by taking the inner product of the signal with each basis sine (i.e. un-normalised correlation). However, that doesn't answer the philosophical question of why don't we just divide by N before returning the output?

Instead, most algorithms divide by N when re-synthesising.

This seems counter-intuitive to me, and (unless I'm missing something) it makes all explanations of the DFT very confusing.

In every scenario I can dream up, the actual magnitude (not the magnitude * N) is the value I need from a DFT operation, and the normalised magnitude is the value I want to input into an IDFT operation.

Why isn't the DFT defined as DFT/N, and the IDFT defined as a simple sum of normalised-magnitude sinusoids?

Whether you scale the output of your DFT, forward or inverse, has nothing to do with convention or what is mathematically convenient. It has everything to do with the input to the DFT. Allow me to show some examples where scaling is either required or not required for both the forward and inverse transform.

Must scale a forward transform by 1/N.

To start with, it ought to be clear that to analyze a simple sine wave, the length of the transform should be irrelevant, mathematically speaking. Suppose N=1024, Freq=100 and your signal is:

f(n) = cos(Freq * 2*Pi * n/N)

If you take a 1024 point DFT of f(n), you will find that bin[100] = 512. But this isn't a meaningful value until you scale it by N. 512/1024 = 1/2 and of course, the other 1/2 is in the negative spectrum in bin[924].

If you double the length of the DFT, N=2048, the output values would be twice those of the 1024 point DFT, which again, makes the results meaningless unless we scale by 1/N. The length of the DFT should not be a factor in these sorts of analysis. So in this example, you must scale scale the DFT by 1/N.

Must not scale a forward transform.

Now suppose you have the impulse response of a 32 tap FIR filter and want to know its frequency response. For convenience, we will assume a low pass filter with a gain of 1. We know that for this filter, the DC component of the DFT must be 1. And it should be clear that this will be the case no matter the size of the DFT because the DC component is simply the sum of the input values (i.e. the sum of the FIR coefficients).

Thus, for this input, the DFT is not scaled by 1/N to get a meaningful answer. This is why you can zero pad an impulse response as much as you want without affecting the outcome of the transform.

What is the fundamental difference between these two examples?

The answer is simple. In the first case, we supplied energy for every input sample. In other words, the sine wave was present for all 1024 samples, so we needed to scale the DFT's output by 1/1024.

In the second example, by definition, we only supplied energy for 1 sample (the impulse at n=0). It took 32 samples for the impulse to work its way through the 32 tap filter, but this delay is irrelevant. Since we supplied energy for 1 sample, we scale the DFT's output by 1. If an impulse were defined with 2 units of energy instead of 1, we would scale the output by 1/2.

Must not scale an inverse transform.

Now let's consider an inverse DFT. As with the forward DFT, we must consider the number of samples we are supplying energy to. Of course, we have to be a bit more careful here because we must fill both the positive and negative frequency bins appropriately. However, if we place an impulse (i.e. a 1) in two appropriate bins, then the resulting output of the inverse DFT will be a cosine wave with an amplitude of 2 no matter how many points we use in the inverse DFT.

Thus, as with the forward DFT, we don't scale the inverse DFT's output if the input is an impulse.

Must scale an inverse transform.

Now consider the case where you know the frequency response of a low pass filter and want to do an inverse DFT to get its impulse response. In this case, since we are supplying energy at all points, we must scale the DFT's output by 1/N to get a meaningful answer. This isn't quite as obvious because the input values will be complex, but if you work through an example, you will see that this is true. If you don't scale by 1/N you will have peak impulse response values on the order of N which can't be the case if the gain is 1.

The four situations I've just detailed are end point examples where it is clear how to scale the DFT's output. However, there is a lot of gray area between the end points. So let’s consider another simple example.

Suppose we have the following signal, with N=1024, Freq=100:

f(n) = 6 * cos(1*Freq * 2*Pi * n/N)  n = 0 - 127
f(n) = 1 * cos(2*Freq * 2*Pi * n/N)  n = 128 - 895
f(n) = 6 * cos(4*Freq * 2*Pi * n/N)  n = 896 - 1023

Notice the amplitude, frequency, and duration differences for the three components. Unfortunately, the DFT of this signal will show all three components at the same power level, even though the 2nd component has 1/36 the power level of the other two.

The fact that all three components are supplying the same amount of energy is obvious, which explains the DFT results, but there is an important point to be made here.

If we know the duration for the various frequency components, then we can scale the various frequency bins accordingly. In this case, we would do this to accurately scale the DFT's output:
bin[100] /= 128; bin[200] /= 768; bin[400] /= 128;

Which brings me to my final point; in general, we have no idea how long a particular frequency component is present at the input to our DFT, so we can't possibly do this sort of scaling. In general however, we do supply energy for every sample point, which is why we should scale the forward DFT by 1/N when analyzing a signal.

To complicate matters, we would almost certainly apply a window to this signal to improve the DFT's spectral resolution. Since the first and third frequency components are at the beginning and end of the signal, they get attenuated by 27 dB while the center component gets attenuated by only 4 dB (Hanning window).

To be clear, the output of the DFT can be a pretty poor representation of the input, scaled or not.

In the case of the inverse DFT, which is usually a pure mathematics problem, as opposed to the analysis of an unknown signal, the input to the DFT is clearly defined, so you know how to scale the output.

When analyzing a signal with a spectrum analyzer, analog or FFT, the problems are similar. You don't know the power of the signal displayed unless you also know its duty cycle. But even then, the windowing, span, sweep rates, filtering, the detector type, and other factors all work to goof the result.

Ultimately, you have to be quite careful when moving between the time and frequency domains. The question you asked regarding scaling is important, so I hope I have made it clear that you must understand the input to the DFT to know how to scale the output. If the input isn't clearly defined, the DFT's output has to be regarded with a great deal of skepticism, whether you scale it or not.


The result if order of two filter are reversed Apply the Composite Laplacian Filter first, then apply the gaussian filter.
Apply the gaussian filter first, then apply the composite laplacian filter.

My work as below:
Here we assume the original image is function $f(x,y)$
$$ (1): (f-\nabla^2f)\ast G = f \ast G - \nabla^2f \ast G$$
$$ (2): (f \ast G) - \nabla^2(f \ast G) = f \ast G - f \ast \nabla^2G $$
I am wondering the statement $\nabla^2f \ast G = f \ast \nabla^2G$ is true or not.

Formally, (linear) derivatives and convolutions may commute, as explained on Wikipedia-Convolution/Properties/Differentiation. This is a major "operation-saving" property: if one wants to differentiate many images and convolve then with a fixed kernel (here a Gaussian), one can convolve the images with the differentiate kernel.
You don't have to differentiate all images, just the kernel once.
Caveat: what is written above could be taken with care, notably when derivatives don't exist theoretically, or with subtleties in implementing derivatives.


script to check zero crossings to check if microphones phases are aligned I have 'n' microphones placed in square fashion and I want to make sure that all the channels must be aligned in time exactly the same if the signal is equidistant from all the 'n' microphones i.e. in the center of the square.

I have written below script to do the difference in zero crossing timings and if the difference is above some precision then print that and bail out.

from scipy.io import wavfile
import numpy as np
import argparse

parser = argparse.ArgumentParser(description='diff the zero crossing of two files')
parser.add_argument('-f1', '--file_name_1', help='provide first file name')
parser.add_argument('-f2', '--file_name_2', help='provide second file name')
parser.add_argument('-p', '--precision', help='precision to compare against', type=float, default=0.0001)

args = parser.parse_args()
print(args)
files = []
files.append(args.file_name_1)
files.append(args.file_name_2)

sampling_rates = []
signals = []
for file in files:
  fs, signal = wavfile.read(file)
  signal = signal / max(abs(signal))                        # scale signal
  sampling_rates.append(fs)
  signals.append(signal)
  assert min(signal) >= -1 and max(signal) <= 1
  print 'fs           ==> ', fs, 'Hz'                       # sampling rate
  print 'len(signal)  ==> ', len(signal), 'samples'

sampsPerMilli = 0
#files should have same sampling rates and length
prev_rate = sampling_rates[0]
for rate in sampling_rates[1:]:
  cur_rate =  rate
  sampsPerMilli = int(rate / 1000)
  if prev_rate != cur_rate:
    print("rates doesn't match %d %d"% (prev_rate, cur_rate))
    exit(0)
  cur_rate = rate

#signal length also should be same
prev_length = len(signals[0])
for signal in signals[1:]:
  cur_length = len(signal)
  if prev_length != cur_length:
    print("length of signals doesn't match for %d %d"% (prev_length, cur_length))
    exit(0)
  cur_length = prev_length

zccs = []
for signal in signals:
  zcc = []
  DC = np.mean(signal)
  newSignal = signal - DC
  for i in range(1, int(len(newSignal))):
    if((newSignal[i] * newSignal[i-1]) < 0):
      #print("crossing at %f seconds"% ((i/sampsPerMilli) * 0.001))
      zcc.append((i/sampsPerMilli) * 0.001)
  zccs.append(zcc)

for a, b in zip(zccs, zccs[1:]):
  if len(a) != len(b):
    print("length doesn't match %d %d"% (len(a), len(b)))
    for c, d in zip(a, b):
      if c-d > args.precision:
        print("precision %f c %f d %f exceeded"% (args.precision, c, d))
        exit(0)

Is there any better approach or can this script be improved?

A better approach for use in the presence of noise and distortion is to use the Wiener-Hopf equations which will provide a least-mean-square solution of the effective "channel" between the microphones. The group delay can be determined using scipy.signal.group_delay.

Further details on this approach specific to microphone captures including Matlab/Octave source code (that can be easily ported to Python Scipy) is on this post:

Compensating Loudspeaker frequency response in an audio signal

In that post, the right channel was equalized to the left channel by passing it through the derived compensation for the effective channel between the microphones, but the post details how you can determine the channel itself rather than the compensation.

You would want to do this with sounding signals that fully occupy the usable spectrum of interest (chirps or white noise generators would be good choices) since the channel can only be determined accurately at the frequencies where signal energy exists. You can do this using one received channel compared to the other as I have done in the link when you don't have a copy of the source, or you can do each left and right channel separately compared to the source if you have the exact waveform that was used to create the source sound.

By using this to derive the channel, and then using the group delay on the channel frequency response, you will see the delay for each frequency in your signal (if there is significant variation, this is group delay variation which is a source of distortion (so you can also use this technique to equalize your room when you do have the source signal). To match the distance between the microphones, you would compare the average of the group delay over frequency, and you can weight this average on frequency ranges that are more important to you.


How to find ROC of system when input is two sided and output is one-sided Given the $\mathcal Z$-transform of input $x[n]$ and output $y[n]$, how can I find the ROC of the system function $H(z) = Y(z)/X(z)$? I have
$$X(z) = \frac{2z\left(z-\frac{10}{3}\right)}{\left(z-\frac 13\right)(z-3)},\quad\text{ROC}: \frac 13 < |z| < 3$$
and
$$Y(z) = \frac{-1}{4}\frac{z}{\left(z-\frac 13\right)\left(z-\frac 14\right)},\quad\text{ROC}: |z| > \frac 13$$

If I substitute $\hat{X}(z) = 1/X(z)$, then $ H(z) = Y(z)\hat{X}(z)$. Then, I know that the ROC of $H(z)$ contains the intersection of the ROCs of $Y(z)$ and $\hat{X}(z)$. But how can I find the ROC of $\hat{X}(z)$? I don't know if the corresponding sequence $\hat{x}[n]$ is causal or not.

The transfer function is

$$H(z)=\frac{Y(z)}{X(z)}=-\frac18\frac{z-3}{\left(z-\frac14\right)(z-\frac{10}{3})}\tag{1}$$

It has three possible ROCs:

$|z|<\frac14$ (left-sided, unstable)
$\frac14<|z|<\frac{10}{3}$ (two-sided, stable)
$|z|>\frac{10}{3}$ (right-sided, unstable)

Only the second ROC, corresponding to a two-sided and stable impulse response, overlaps with the ROC of $X(z)$. The intersection, and, consequently, the ROC of $Y(z)$ would be $\frac13<|z|<3$ (equal to the ROC of $X(z)$). However, the pole at $z=3$ is removed by pole-zero cancellation, which results in the final ROC $|z|>\frac13$ for $Y(z)$.


Precisions on Hilbert Huang transform After investigating signal analysis methods based on empirical mode decomposition (EMD), I found that recent developments are mostly related to the Hilbert Huang transform (HHT) and the Local Mean Decomposition (LMD) method.

I have been reading a few articles on the subject and I would like to have your opinion regarding HHT.

EMD leads to intrinsic mode functions that are particularly well suited for Hilbert transform. HHT seems to be really widely used for many distinct types of industrial or academic applications.

Am I right to think that HHT is essentially EMD + Hilbert transform? Do you consider the specificity of the HHT lies in the sifting process based on spline interpolation (to the contrary of the LMD using a moving average algorithm for instance)?

After going through the literature regarding HHT and EMD, I found that the "Huang" part of HHT comes from the fact that he is the one who proposed EMD in the first place. That explains the name of the method...

For more developments regarding EMD and HHT, I recommend the papers by Rilling et al. "On empirical mode decomposition and its algorithms". For the lucky ones who speak French, Rilling's PhD. thesis on EMD no longer seems to be available online; the document seems very complete and contains a very detailed mathematical analysis of EMD. Related articles are also available in English here on Google Scholar.

One may sum-up HHT this way:

EMD: decomposes the initial signal as a list of Intrinsic Mode Functions (IMF);
Hilbert tranform computes the instantaneous frequencies associated to the IMF (which are precisely well suited for such transformation)
Hilbert spectrum, meaning a representation of the amplitude of the IMF in a frequency/time domain using the instantaneous frequencies.

Here is a simple example based on the time signal: $y(t)=t^2sin(t^2)$ for $t\in[0;~16]$s. This signal features a quadratic growth of its amplitude combined with a linear growth of its frequency. Here is what a HHT of the signal leads to:

signal of interest

IMF obtained from EMD

$y(t)$">

instantaneous frequencies

Hilbert spectrum (white to black from 0 to maximum amplitude)

Some peaople see the HHT as a generalized Fourier Transform in the sense that the decomposition of the signal interest by HHT leads to both amplitude and frequency varying time signals.

A signifiant drawback of HHT lies in its sensitivity to edge effects (what happens for the signal close to its left and right boundaries in 1D). Several techniques exist in order to mitigate these effects. Articles by Rilling mentioned above go with mirroring strategies while other engineering oriented strategies may involve waveform matching techniques.


Beginner level : Help with terminologies : smooting length, length of channel, channel equalization delay, blind system identification, equalization Say the channel model is an univariate FIR filter with true coefficients $h=[1,a_1, a_2 ]$. I am learning algorithms for system identification and channel equalization. For this, I am implementing programs provided in Matlab. One such Algorithm is Constant Modulus Algorithm for blind channel equalization.

Code is the link for the code : lines 80-82 present the variables $L,ChL,EqD$  which denotes the smoothing length, length of the channel & channel equalization delay respectively. I have confusion regarding the terms. Shall appreciate explanation with an example pertaining to my case.

1. Smoothing length : In text books, $L$ stands for the channel order, which for my model, Moving Average, MA(3), the order then = 3 ( I guess including the 1). FOr Autoregressive model AR(p) $p$ is the model order. I am confused if smoothing length has got anything to do with Kalman Smoother ? What is the difference, if any, between smoothing length and channel order? What will $L$ be for the FIR model that I have?

2. Length of the Channel, ChL = Is this related to number of data samples ?

3. Equalization delay, EqD = In general, I was under the impression that the model order is determined by the number of delays. So for MA(q) : $q$ is the number of delays, similarly for AR(p) : $p$ is the number of delays. Then what is Equalization delay?

When considering blind algorithm, we assume no knowledge of the process noise that drives the model. But, When we assume that the model order is known, then does it still become blind or non-blind estimation / equalization technique?

Please correct me where wrong.

The term smoothing length (also called observation interval) is simply the length of the equalizer, i.e. its number of filter taps. In the code the filter length is L+1, so L is the filter order. Similarly, the length of the channel refers to the length of the FIR filter simulating the dispersive channel. Again, the channel filter length is actually ChL+1. What is called equalization delay in the code is the total delay introduced by the channel and by the equalizer. It is needed for aligning transmitted and received (and equalized) symbols.


Removing known sequence from an audio sample I have an audio recording of a conversation which is obscured by a movie playing on a television.

I have an audio recording of that same movie playing on that same television with no conversation in the background.  The recording device is identical as is its placement in relation to the speakers.  The volume is close to identical.

How do I remove the audio from the television so that I can unobscure the conversation?

Obviously ambient noise is different and system/thermal noise in the recording device will be different.  The phase of the two audio recordings will be different as there is some time error bounded by $1/F_s$.  ($F_s = 44100\textrm{ Hz}$)

I was unsuccessful in attempting to find a time offset between the two signals to perfectly align the audio.  I loaded the files in MATLAB and took the sum of the absolute difference between the two samples for multiple offsets and then selected the one with the lowest error.

The reconstructed signal being the difference between the shifted reference audio and the conversation.

 reconstructed_signal = conversation - reference_audio(offset:(length(conversation)-offset+1));

This time domain approach has not worked.  I am considering that I should do a difference of the frequency data of the two samples with a sliding window.

Do you think that this would be a good approach?

From the Description, I think Adaptive filters will work out for this application.
As you have the reference signal("I have an audio recording of that same movie playing on that same television with no conversation in the background") and the signal from which you need to remove this reference.


Real-time oscillation-like noise filtering I am making real-time measurements, which are affected by some noise. This noise has oscillating pattern, i.e. the total signal tends to oscillate around the mean value. The examples of measurements are shown in figures below.

Filtering this noise after the signal is received is quite straitforward - we may use mean filter, for example.

What I would like to do is to implement real-time noise filtering, so that filtered signal would be close to its constant mean value.

What kind of methods of filtering could you recommend? I thought about Kalman filter, but the state model is unknown. Maybe there is something else?..

P.S. I am totally unexperienced in signal processing, please excuse me for any ignorance.

UPD: I am not considering mean filters, because the two pictures below represent only particular stages of evolution. It means that the signal evolves in an unpredictable manner in between. Mean filter usually leads to loss of information about such evolution.

Marcus and Ben have very clearly identified the problem: the signal is simply too small relative to the least significant bits (LSBs) of the "virtual machine"'s ADC. The sensible solution is to increase the signal by use of gain. From the provided plots, a signal gain around 100 would be a good starting value. Alternatively, the ADC resolution would have to be increased (> 4 more bits), so the LSBs would be smaller. In a real system, the first option would be the only feasible option, since amplifiers are relatively inexpensive and high resolution ADCs are most certainly not. Since the goal is real-time noise filtering, rather than Lazarus filtering, i.e., post-mortem reversal of fatal mis-quantization, dithering and signal averaging is ruled out.


Sharp EL-5500III (PC-1403) cassette tape encoding I have an old audio cassette tape of digital data I converted to .wav and am hacking to recover the original data. I demodulated the 4K/2K FSK modulation with minimodem. Then I get stuck with insufficient understanding of the NRZI bit-stuffing algorithm. Short-cut: does anybody know of sox-like or minimodem-like program that does NRZI decoding through Linux pipes?

Many references explain that bits are stuffed into NRZI data streams after 5 consecutive bits of one polarity. But to be precise, I haven't found a clarification. Does this mean:

a bit stuff happens after 5 bits, or
a bit stuff happens in order to prevent a run of 6 bits?

By example, if the original data is ...11111001... does it get transmitted as is, or does a bit stuff occur to send the data ...111110001... ? If the bit-stuff is only to cause a transition after 5 bits, then the stuff is not necessary. But the docs say "after 5 bits" which would create unnecessary stuffs if the original data had inherent runs of exactly 5 bits.

Here's the bit stream I'm trying to decode.

root@kali:~# minimodem --rx -R 44100 --mark 4000 --space 2000 --binary-raw 26 -f monoshort8008.wav 500
CARRIER 500 @ 4000.0 Hz
11111111111111111111111111
11111111111111111111111111
11111111111111111110111010
00001111100000100000111110
00001000001111100000100000
11111000001000001111100000
10000011111000001000001111
10000010000011111010101011
11111110100010001011111000
00100000110000010010111000
00101100110101110111111000
10101000110000010101111000
00100000110100010001011000
00100101110101110101011000
10101001110110010101111000
10101000110101110000011011
00100000110101110100011001
00101011110110010100011000
00101011110000010000011010
00100111110000010110011010
11101011110001010100111000
00101011110000010000011001
00100001110000010110011010
11101001110001010100111000
00101011110000010000011011
00100100110000010110011000
11100010110001010100011000
00101011110111110111111011
11101111110000110111011110
NOCARRIER ndata=29 confidence=150.653 ampl=0.126 bps=500.04 (0.0% fast)

Notice the sets of 5 bits near the beginning. Notice also the periodicity of 13 bytes, visible as vertical column runs of 1 or 0 through all the data.

Andy suggested Octave check the frequencies of the modulation.  Here is a plot of the data - x axis is an unscaled count of data points. Data was collected at 44,100 Hz.  Working the scale factor shows 2KHz and 4KHz, 2msec bit time. Matches Audacity waveform view and Minimodem FSK demodulation.


packet encoder RF packets by using LPF for determining if it's zero or one I'm transmitting RF packets with frequency 868MHZ ,
symbol rate (it's the same as bit rate) is 50KBaud , and Im sampling them by a dongle (fs=20480KHZ) to capture those RF packets, the captured samples are IQ RAW DATA and I save them into a file , by matlab Im using this file and plot all the captured data(packets), I see the RF packets as peaks where there's packet found, otherwise noise.
I have built demodulator which the output of its values ranges between 1 and -1 , so I have arrived to a file that has values that ranges between 1 and -1 , this demodulated IQ samples I want to convert them to binary 0 or 1 , so Im trying to build encoder in Matlab code by using LPF , how could I do that? for more clarification I want to input to my LPF the demodulated IQ samples and then the output would be something like:

I want this because if I have the output of LPF then I can walk through that output and decide where I have "zero-Crossing" this means when there's a step (zero-crossing) then the value would be changed from zero to 1 or 1 to zero (depend on the value what was before) ,
how can I do that? the cutoff frequency of my LPF  must be my bit rate frequency right? how can I also do LPF in matlab? thanks alot !!!!
second question, lets assume I have the output of my LPF , then how can I walk through the output in order to decide zero/one (deciding zero crossing-high to low or low to high-), this part called FLOW ESTIMATION , how can I do that?

If the constellation at the tarsmitter is symmetric about 0 on the real line then you can simply threshold the output between -1 and +1 with 0 as the threshold. So anything above zero is a logical 1 and everything below zero is a logical 0.
This detection is rule is optimal when the noise is additive.


Digital frequency for images From this Wikipedia link, the fastest fequency for digital frequency is 0.5cycle/sample. It can be rewritten as 1cycle/2sample.

What does it mean?
Does it mean that at least two pixles are required to form one cycle?
Is there any requirement on the values of two samples?
One must be 0 and the other must be 255? If so, why?

Edit:
The reason I have such question is that different cutoff frequencies result in different blurring degree for filtering in frequency domain.

From the link:http://www.imagemagick.org/Usage/fourier/#fft_partial

The figure shows the comparison. I hope to estimate blurring degree from cutoff frequency in digital frequency. Maybe I can estimate corresponding standard deviation of Gaussian filter from cutoff frequency.

Highest digital spatial frequency is 0.5 cycles/sample. It's related to Nyquist criterion.

  Does it mean that at least two pixels are required to form one cycle?

Yes. Minimum two pixels are required to form the cycle.Try to fit a one cycle of sine wave in minimum number of pixels. You will find that two minimum pixels are needed to fit the sine wave without the loss. Anything less than two pixel you can see the loss.Finally if you fit a full sine wave in single pixel,the pixel is average of whole cycle resulting in total loss.

  Is there any requirement on the values of two samples?

No, there are no requirements like that.

For the edit

The blurred image after Gaussian smoothing is given by $g(x,y) = f(x,y)*h(x,y)$
,where f(x,y) is the input image, h(x,y) is the gaussian smoothing.

In frequency domain G(u,v) = F(u,v)H(u,v)

H(u,v) = G(u,v)/F(u,v)

The Gaussian blurring function h(x,y) = ifft(H(u,v)) .

Using only the cutoff frequency from the FFT of the blurred image g(x,y) you will get the good approximate of your gaussian filter but not accurate one. The reason behind this is the input itself might not contain high frequencies but from the FFT you might think of it caused by the Gaussian smoothing. But for most natural scenarios you can get good approximate of your Gaussian smoothing from cut-off frequency itself.


FIR Filter design by non-equal samples I'd like to design a Digital Differentiator with a pre-specified lowpass characteristic.

I've tried Wiener Filters (using inverseFFT to transform from frequency to time-domain), Savitzky Golay and others. Intuitively, I've come to the conclusion, that by specifying DFT points, i.e. uniformly distributed samples of the frequency response, one can achieve only a maximal low-pass frequency. For example using $N$ points, the minimal lowpass frequency would be $1/N$ bins of the FFT, thus $\pi/N$ $rad/sample$.

Is this notion correct?

And further along this path: if it's true, is this a universal boundary, or are there methods, where one can implicitly or explicitly specify not-equally sampled frequency points and thus achieve lower stopband frequencies with FIR filters?

For IIR filters I think this is no problem, since with a simple $1/(z+a)$ one can choose $a$ to achieve a variable stopband frequency just as in the continuous case.

The Remez exchange method isn't limited to DFT bin constraints so to answer your question, no FIR filters do not have the limitation you hypothesized.

If you can find METEOR, a program written in part by Tom Parks, you can be fairly free in your filter specifications provided they are feesable

and there is

http://gmeteor.sourceforge.net/


Efficient algorithm for computing area under convolution? Consider $N$ discrete signals $x_1(n), x_2(n), \ldots,x_N(n)$ each a bounded support of size $M$.

To convolve them, we can zero-pad each of them, multiply their FFTs, and take the result's inverse FFT.
I believe this requires taking $N$ FFTs of length $N(M-1)+1$ each, taking $O(MN\,\log(MN))$ time.

However, what I need isn't the full convolution. I only need the area under the convolution above zero:
$$\sum_{k=0}^\infty y(n)$$
assuming $y$ is the convolution of our signals:
$$y = x_1 * x_2 * \ldots * x_N$$

Is there any algorithm that allows me to compute this value with a better time complexity than that required for computing the full convolution, or are the problems equally difficult?

Note that

$$\sum_ny[n]=Y[0]$$

where $Y[k]$ is the DFT of $y[n]$. Since

$$Y[k]=\prod_iX_i[k]$$

with $X_i[k]$ the DFT of $x_i[n]$ (appropriately zero-padded), you have

$$\sum_ny[n]=Y[0]=\prod_iX_i[0]$$

And since

$$X_i[0]=\sum_nx_i[n]$$ you get the final result

$$\sum_ny[n]=\prod_i\sum_nx_i[n]$$

the computation of which requires $N-1$ multiplications and $N(M-1)$ additions.

This little Matlab script illustrates how it works:

M=100;          % length of each signal
N=10;           % number of different signals
X=randn(M,N);   % each column of X contains a signal
Xs=sum(X);      % sums over each signal
Xsp=prod(Xs);   % product of the sums = desired result

% compute convolution of N signals
y=X(:,1);
for i=2:N, y=conv(y,X(:,i)); end
ys=sum(y);      % this should equal Xsp

[Xsp, ys]       % .. and it does


Solving for impulse response I'm working on some problems that instruct me to find the impulse response of a given system. For example

$$y(n)=4y(n-1)-4y(n-2)+x(n)-x(n-1)$$

The first step I take is to find the homogeneous solution of the equation like so:

$$ \lambda^n-4\lambda^{n-1}+4\lambda^{n-2}=0
\\\lambda^{n-2}(\lambda^2-4\lambda+4)=0
\\\lambda^{n-2}(\lambda-2)(\lambda-2)=0
\\\lambda = 2,2
\\y_{h}(n) = C_{1}(2)^n + C_{2}(2)^n
 $$

Is this correct homogeneous form for when you have two of the same values for lambda? I know it is when the values are different from each other.

I also have another question about this practice. When I've finished this step, I then plug n=0 and n=1 into the original equation to find the values of C1 and C2. How do I know how any times to plug a value into the original equation?

Thanks.

Your characteristic polynomial is indeed

$$\lambda^2-4\lambda+4=(\lambda-2)^2$$

If you had two distinct roots $r_1$ and $r_2$ the general solution would have the form you suggested:

$$y(n)=C_1r_1^n+C_2r_2^n$$

However, in your case you have a double root $r=2$ which results in a general solution of the following form:

$$y(n)=C_12^n+C_2n2^n\tag{1}$$

The constants $C_1$ and $C_2$ are chosen to match the initial conditions and the input signal $x(n)$, which in order to compute the impulse response is a delta impulse $\delta(n)$. Assuming $y(-1)=y(-2)=0$ we get

$$y(0)=x(0)-x(-1)=1\\
y(1)=4y(0)+x(1)-x(0)=4+0-1=3$$

Comparing to (1) gives

$$y(0)=C_1=1\\
y(1)=2C_1+2C_2=3\Rightarrow C_2=1/2$$


RMS values above 0 dbFS? I wonder why several RMS meters I'm using for music production show RMS values above 0.

For example, when using an input signal with a constant DC of 1 (0 dbFS), I get these results, although Wikipedia says that in such cases, the RMS is equal to the amplitude of the direct current.

Samplitude Pro X (RMS: 2.9)

TT Dynamic Range Meter (RMS 3.0)

Whereas other meters show the values I expect:

Voxengo SPAN (RMS: 0.0)

Why is it that meters come to such different results?

0 dBFS refers to a full-scale sine wave or more rarely a full-scale square wave.
ITU-T telephony Recs. P.381 and P.382 and reportedly the paywalled digital audio device measurement standards AES17-1998 and IEC 61606 define that full-scale sine wave is 0 dBFS. ITU-T Rec. P.10/G.100 (11/2017) does not make the distinction and defines dBfs (sic) as "relative power level expressed in decibels, referred to the maximum possible digital level (full scale)".
Full scale DC has the same RMS value as a full-scale square wave, but is 3.010299956 dBFS RMS if the 0 dBFS reference is a full-scale sine wave. Looks like the different programs use a different reference (square or sine). Samplitude Pro X might use some approximation of RMS that gives in this case a 0.1 dB unit error.
The calculations:
$$\int_0^{2\pi}\frac{\sin(x)^2dx}{2\pi} = \frac{1}{2}$$
$$10\log_{10}\left(\frac{1}{2}\right) \approx -3.010299956$$


Question on PCM Sampling and Quantization order I had a doubt in my mind regarding the orders of sampling and quantization in PCM. What is the impact if the order is reversed, that is, the continuous time signal is first quantized and then sampled? I find a source on the internet claiming that quantization noise will increase with the reversal (link given), but I'm not sure why. Personally, I feel it should be unaffected as we're simply quantizing as in allowing a finite set of amplitudes first and then sampling those amplitudes. Am I missing something?
Link which claims the same: (Question 16)
https://electronicspost.com/signal-and-systems-interview-questions-and-answers/

That link provides a very superficial Q & A on signals and systems.
It's not clear which kind of a sampler or quantizer they are referring to.
However, if sampling and quantization are defined as ideal mathematical operators then,
you should get exactly the same result by altering their orders. Hence the sample values should be same (meaning that the same noise also).
May be they're referring to some unstated practical ADC schemes (that I cannot recall) which do not perform ideal operations and might create different results.


Discrete Fourier Transform bandwidth? In a problem I'm working on, I have a signal that is sampled 256 times for 3.5 microseconds.

After performing a Discrete Fourier Transform on the data, what is the frequency step between the calculated values of Cn? The prompt states that this is also called the bandwidth of each point in the power spectrum.

How can I go about solving this problem?

Thanks

Your specification of the particulars for your situation aren't clear to me, but I'll explain how you can arrive at the answer. The DFT is defined as:

$$
X[k] = \sum_{k=0}^{N-1}x[n] e^{-j2 \pi nk/N}
$$

The (normalized) frequency associated with the $k$-th bin is $\frac{2\pi k}{N}$. Therefore, the spacing between the $k$-th and $(k+1)$-th bins is $\frac{2\pi}{N}$. Recall that the conversion between unnormalized frequencies (i.e. before sampling) and the normalized angular frequency $\omega$ used in the DTFT is:

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

That is, a time-referenced frequency $f$ before sampling maps to $\omega$ in the normalized frequency domain, where $f_s$ is the sample rate of the signal (i.e. how many samples are taken of the signal per second). You can then map the bin spacing above into unnormalized frequency:

$$
\Delta f = \frac{2\pi}{N} \frac{f_s}{2\pi} = \frac{f_s}{N}
$$

This simple relationship just says that the spacing between bins is equal to your signal's sample rate divided by the DFT length. It's unclear from your question whether the interval between each sample is 3.5 microseconds or whether your block of 256 samples covers an overall duration of 3.5 microseconds. Based on whichever is true, you can calculate the appropriate value.


How can I calculate the entropy of a signal that's not independent from itself? I really don't know how to word this better, sorry if that title is confusing.

Basically, I'm wanting to use the noise from an AM radio audio signal to generate random numbers. To do this effectively, I need to know how many bits of entropy I can expect from some number of bits of signal data.

If every sample were independent, I could easily calculate the entropy using probabilities. The problem is that the samples are dependent on the samples that came before it and the samples that come after. I don't know how to account for this.

How can I calculate the entropy (in bits) of a section of my signal?

I think that the answer is that there is no correct answer.  Whilst this is a simple question, the answer is kinda open.

You're asking about the entropy of correlated data with respect to random number generation.  The key reference for this is NIST Special Publication 800-90B, Recommendation for the Entropy Sources Used for Random Bit Generation. Look to sections 5 and 6. You be faced with two questions; Is my data correlated (non IID) and if so, what is it's entropy rate? You have already realised that the standard $ -K\sum_{i=1}^n {p_i \log p_i} $ isn't appropriate.

Is your data correlated?

This is the first step, and NIST suggests permutation testing, as below:-

There are 11 tests NIST recommends, all of which are somewhat complex.  To help (not) they have release some code on GitHub that implement these tests.  I'm unsure of how official this is, but see the well received How to interpret the entropy results for a NIST test file as to their soundness. An immediately apparent problem for implementation is that as $ T, T_i \in \mathbb{R} $, the equality test in Fig.4, section 2.2.2 becomes difficult without bounds. In a nutshell, the code doesn't work well enough for field use.

What is the entropy?

If you manage to satisfy yourself that the data is non IID, NIST recommends taking the lowest value of another 10 complex statistical tests as the min. entropy value.  Similarly, the Python and C++ code  has problems outputting a coherent and credible entropy estimate. Perhaps you can implement the tests yourself in another fashion.

What I do.

I too build random number generators and I use similar permutation tests, but focused solely on a compression test. Compression exploits correlation and eliminates redundancy.  Compression is also used in both the NIST non IID confirmation and entropy estimation tests above.

One can calculate what I have termed a correlation factor, $ c_f $, explored in What exactly does compression say about correlation of data?  If $ c_f \simeq 1 $ then the data is IID, and non IID if  $ c_f > 1 $. I then use the non permuted compressed file size and divide twice by two to get a conservative entropy rate.  This method seems to be as good as any, and more robust that most. I estimate about 1.6 bits / byte from a 24V Zener diode at 10kSa/s. That's 16kbps of cryptographic strength entropy.

Notes.

I now use fp8 and paq8px compressors as they are about the best I can find that will compile for me. They'll compress IID data to within 1% of the theoretical Shannon limit.
Divide by two as a security precaution to allow for improvements in compression algorithms, whilst realising that these improvements are asymptotic with diminishing returns.
Divide by two again just because.


Fast Morphological Erosion / Dilation with Variable Window Size I want to perform the following operation on an image: consider a square window around every pixel, find the maximum value in this window, and replace the central pixel value with this maximum value. If I consider a $n \times n$ square window around every pixel, this becomes the morphological dilation operation (with a square $n \times n$ structuring element).

In my work I need to have windows of different sizes around every pixel, i.e. I might need a $n \times n$ window around pixel $A$ and a $m \times m$ window around pixel $B$.

Questions:

(1) Does this operation have a name? Is it widely used?

(2) Do you know of any fast algorithm to perform this operation?

Relevant info: For fixed window sizes, there does exist a fast algorithm, for example Algorithm 1 in the following paper (which in fact requires just a constant number of operations per pixel, regardless of window size):

Chaudhury, K. N. "Acceleration of the Shiftable O(1) Algorithm for Bilateral Filtering and Nonlocal Means." IEEE Transactions on image processing, 22.4 (2013): 1291-1300.

The trick when using masks is that if the window is shifted one to the left, only one row is added and one removed. Therefore one can e.g. calculate the mean of a 5x5 window quicker, as if all 25 values needed to be added and divided again.

In erosion or dilation this fact could be used to, if one would store the max of each row. Then only the maxes of one new row would be needed and a comparison between all rows. However I guess nobody does that, since the overhead is quite big and big windows would need additional space for each row's max...

Your problem however couldn't even use this, since the window size is changing. I assume the changes are random. Therefore there will be no quicker way than moving from pixel to pixel and apply the different window sizes.

But why would someone use different windows for different pixels?


How to recover a signal from its complex envelope I know there's similar questions here, but I did not fully understand the answers I have read:
in MATLAB I consider the complex envelope of an output signal:

y_low_out   = A_low_out .* exp(1i*Phi_low_in) .* exp(1i*Phi_low_out);

Now, I want to recover a sinusoidal signal, what does y_out look like in terms of superpositions of sine-waves? My guess is:

y_out = A_low_out .* sin(2*pi*f_c*t + Phi_low_in + Phi_low_out)

But this is pure guessing. $f_c$ is the carrier frequency.

Some background to what I am trying to do:
I consider an input signal:

x(t) = A_in .* sin(2*pi*(f_c+f1)*t + Theta1) + A_in .* sin(2*pi*(f_c+f2)*t + Theta2))

This input signal is converted to an equivalent lowpass signal (complex envelope) with amplitude A_low_in and phase Phi_low_in and sent through a nonlinearity which changes both amplitude and phase. The amplitude of my output signal (still complex envelope) is A_low_out, the phase change is Phi_low_out. This description refers to the code that I have written. I have determined A_low_in by taking the absolute value of the analytic signal of x. For some reason A_in and A_low_in are not the same. I figured A_low_in would be twice A_in, but it's not exactly (maybe numerical differences).
Anyways, in the end I get the complex envelope y_low_out. And I want to convert it to y_out.

According to the usual definition, if you have a complex low pass signal (complex enevelope) $u(t)$ and a reference frequency $f_c$, then the corresponding pass band signal is

$$u_p(t)=\text{Re}\{u(t)e^{j2\pi f_ct}\}\tag{1}$$

So if $u(t)$ is given by

$$u(t)=A(t)e^{j\phi(t)}\tag{2}$$

the real-valued pass band signal would be

$$u_p(t)=\text{Re}\{A(t)e^{j\phi(t)}e^{j2\pi f_ct}\}=A(t)\cos(2\pi f_ct+\phi(t))\tag{3}$$


Calculating phase with noninteger periods of a sinusoid How do we calculate the phase of a signal with a noninteger period of a signal? For example, if a signal has a frequency of 1kHz and is sampled at 8kHz, 8 samples (1 period) can be used in the FFT or Goertzel to calculate the phase. However, if a signal has a frequency of 900Hz or 1.1kHz and is sampled at 8kHz, then 8 samples (1.11 or 0.909 periods respectively) will produce incorrect phase calculation. How do we compensate for a noninteger sample period?

Do an fftshift (rotate the input vector by N/2) before doing the FFT.  Or for a phase measurement at a fixed point in time, start sampling the data earlier so the N/2 time sample in the FFT input vendor is at the time point where you want the measurement, and then flip the sign of the phase of every other bin of the FFT result (the odd numbered ones).  This will move the phase reference point to the middle of the data window, where an input sinusoid will not have a discontinuity in phase that changes with frequency.  This also allows interpolation of phase between FFT result bins.

The evenness/oddness ratio, and thus the FFT phase, does not change with frequency when a sinusoid's phase is referenced to (for example both generated at and measured with respect to) the center of the FFT aperture.


Convert a Park McClellan FIR Solution to Achieve Stop-band Roll-off Often it is more desirable to have a roll-off versus frequency of the rejection in the stopband; for example resampling filters where all the alias image locations fold into the first Nyquist Zone. What is a VERY simple change you can make to the coefficient solution provided by the Parks-McClellan algorithm that will convert the "equiripple" stopband to one that rolls off with greater rejection (at the expense of exceeding target limits near the transition, and increasing passband ripple)? Note: This is a great "trick" that I learned from fred harris.

Note- for full credit please explain why this "trick" works.

To help illustrate the end result, see the plot below showing a standard design using P-McC before and after the simple modification.  Here we see the stop band rejection with a roll-off instead of being flat, but increased passband ripple and degraded rejection immediately near the beginning of the stopband, but overall significantly better rejection.

Before (blue) and after (red) modifcation to filter coeff. Note the region close to the transition where the stop band is worst, but overall rejection is much better:

Zoom in showing trade-off of degraded ripple in passband.

That trick must have to do with manipulating the first and the last tap of the Parks-McClellan (Chebyshev) solution. This is also explained in the last paragraph of this answer, which deals with a very extreme case. For equi-ripple designs the first and the last tap tend to be a bit larger in magnitude than what would be expected if one just extrapolated the impulse response. These two extremal taps correspond to time-domain echos and generate the equi-ripple behavior in the frequency domain. By decreasing the size of these two taps, the stop band ripple can be reduced. Any heuristic can be used to do so. I don't know what Harris suggested, but by simply choosing $h[0]=h[1]/2$ and $h[N-1]=h[N-2]/2$ (where $N$ is the filter length) I get an even faster roll-off (albeit a less regular one), at apparently no extra cost (compared to the solution in the question):

Note that this trick only works well for filters with a relatively narrow pass band. Equi-ripple filters with wider pass bands usually do not show that irregularity in their impulse response.


Role of Kalman filter in nonlinear dynamics I am quite interested to know the significance of using kalman filter,unscented kalman filter and extended kalman filter in chaos synchronization when infact the very basics of chaos synchronization(physics perspective) never mentions about using kalman filter!However, these papers KF1 and alike use kalman filter for synchronization. There are various methods of chaos synchronization viz generalized, phase, projective, adaptive Boccaletti_Kurths et al and wikipedia considered to be the bible in thi field, do not mention about using kalman filtering and alike methods. Would be obliged if the role of this technique is explained since I think I am missing on important theory. Thank you

Kalman filters are useful if one is dealing with noisy systems, when one either has no knowledge about the source of the noise or the noise is created by a high-dimensional nonlinear system with many positive Lyapunov exponents.

Chaos synchronization is dealing with small low-dimensional systems which usually posses only one or two positive Lyapunov exponents. Applying a Kalman filter does not make much sense, since many better suited methods exist to deal with such systems, as one knows, that the system is not a stochastic one.


Same z transformed function, but different answers of inverse z transform? Given a $\mathcal{Z}$ transformed function $E(z)=\frac{1}{z+4}$.

I know there are several ways to get the inverse $\mathcal{Z}$ transform of this function :

Using partial fraction

$$E(z)=\frac{1}{z+4}$$

$$\frac{E(z)}{z}=\frac{1}{z(z+4)}$$

apply partial fraction here,

$$\frac{E(z)}{z}= \tfrac14 \left(\frac{1}{z} - \frac{1}{z+4} \right)$$

so $E(z)$ is,

$$E(z)= \tfrac14 \left(1 - \frac{z}{z+4} \right)$$

as you know, it is easy to use inverse Z-transform here.

$$ e[n] = \tfrac14 ( \delta[n] - (-4)^n ) $$

inversion formula method

$E(z)$ has simple pole at $z=-4$, the residue is evaluated as

$$\begin{split}
\text{residue}_{z=-4} &= ((z+4)\times E(z)\times z^{k-1})_{z=-4}\\
&= (1\times z^{k-1})_{z=-4}\\ &= (-4)^{k-1}\\ &= -(1/4)\times (-4)^k \end{split}$$

The problem is that inverse $\mathcal{Z}$ transform of the same function $E(z)$has 2 different answers. As you see, there is no delta function when I used inversion formula method.

At first, I thought it was because a simple pole, $z=-4$, is outside of the unit circle, but same things happen when I used a simple pole which is inside of the unit circle.

Why does this difference happen? What is the real answer?

Here's a much simpler method. We have

$$E (z) := \frac{1}{z+4} = z^{-1} \left(\frac{1}{1 - (-4) \, z^{-1}}\right)$$

which we assume to be the $\mathcal Z$ transform of a causal discrete-time signal. Consulting a table of $\mathcal Z$ transforms, we conclude that

$$\mathcal Z^{-1} \left\{ \frac{1}{1 - (-4) \, z^{-1}} \right\} = (-4)^n \, u (n)$$

where $u$ is the discrete-time step function. Hence, the inverse $\mathcal Z$ transform of $E$ is simply the discrete-time signal above after a unit delay, i.e.,

$$e (n) = (-4)^{n-1} \, u (n-1)$$

which is zero at $n=0$. Hence, your 1st candidate answer is (almost) correct. Your 2nd candidate answer is only correct if $n \geq 1$. When working with functions, do not forget what their domains are.


Get minimum phase from function Why is it that reflecting any poles or zeros of a rational function across the unit circle gives a minimum phase system? Here's an example, it seems reflecting any poles or zeros would result in the zeros being outside the circle?

Why is it that reflecting any poles or zeros of a rational function across the unit circle gives a minimum phase system?

It doesn't. You are starting with a wrong assumption.

Here's an example, it seems reflecting any poles or zeros would result in the zeros being outside the circle?

Correct. Hence it's nonsense and you shouldn't do it. Assuming that the 'x' marks the poles, the system is unstable or non-causal anyway.
Roughly speaking:

A system is minimum phase if all zeros and poles are inside the unit circle.
A system with zeros outside the unit circle can be made minimum phase by reflecting the zeros at the unit circle
If have poles outside the unit circle, chances are something went wrong already and minimum phase is the least of your problems.


Generating data with Gaussian distribution We want to generate a Gaussian vector, (call it $\Delta$) - with an arbitrary size-, with zero mean and variance of $\alpha$. If $\alpha$ is chosen such that $\parallel\Delta\parallel_2 \leq 0.5$ with probability of e.g. $90\%$.

How can we do that in Matlab?

Let $\Delta = [x_1, \dots, x_N]$ be your desired vector where $x_i$ are iid $\mathcal{N}(0,\sigma^2)$. Do you expect that $\mathrm{Pr} \left\lbrace ||\Delta||_2^2 < 0.5\right\rbrace = 0.9$ ?

By definition $y = ||\Delta||_2^2 / \sigma^2 = \sum_{i=1}^N \frac{x_i^2}{\sigma^2} \sim \chi^2_N(y)$ is Chi squared distribution of degree of freedom $N$.

Thus $z = ||\Delta||_2^2 = \sigma^2y \sim \frac{1}{\sigma^2} \chi^2_N(\frac{z}{\sigma^2})$

Denote the CDF $\mathrm{F}(z)$, you want something like $\mathrm{F}(z = 0.5) = 0.9$ with $\sigma^2$ is unknown variable.

A closed form expression of $\mathrm{F}(z)$ is available at wikipedia Chi squared distribution.

In MATLAB you can use the two functions chi2pdf() and chi2cdf() for numerical brute force solution for $\sigma^2$.

And after having the value of $\sigma^2$, just use the function randn() to generate your vector $\Delta$. For example with size N and $\sigma^2$ = sigma2, delta = randn(1,N) * sqrt(sigma2).


What's a "delay line"? I know what a delay is, but
what's a delay line? What's the reason for calling it a "delay line"?

The term delay line comes from (analog) telecommunications where timing between different signals needed to be adjusted by delaying one of them. For example, the chrominance and luminance channels in a cathode-ray-tube colour television.

The reason it's called a line is because that's all it really is: a wire. The website linked to above has the following diagram of one.

In the discrete (digital) domain, I've more usually heard delay lines called tapped delay lines.  That's because each sample is a delay, and to form an FIR filter, we need access (to tap into) the signal at different delays.  The line of $z^{-1}$ boxes in the diagram below is what I would call a tapped delay line.


Design a FIR filter with a Kaiser window Here is the exercise I am working on.

  Use the window design method with the Kaiser window to design a linear phase high pass FIR filter of minimum order that meets the following specification:
  $\omega_P=0.8\pi$,$\omega_S=0.7\pi$, maximum pass band error=0.1dB, minimum stopband attenuation = 75dB ( you can use Matlab but you need to ensure that the final filter meets the spec. Recall that for window design methods, the errors in pass and stopband will be equal. Be sure to take it into account in your answer)

I know how to design a FIR filter using window method. What I do not understand is the last two sentences. What is "the errors in pass and stopband"? What it mean by "take it into account"?

The Kaiser window has a parameter, usually called $\beta$, which is determined by the specified peak error $\delta$. For your specification, you need to determine which of the two requirements (maximum pass band error or minimum stop band attenuation) is more stringent, and choose the allowable peak error $\delta$, and hence the parameter $\beta$, accordingly.


Anyone know what algorithm the Spice AC Noise Analysis uses? Anyone know what algorithm the Spice AC Noise Analysis uses?

http://vision.lakeheadu.ca/eng4136/spice/noise_analysis.html

Is it some spectral modeling synthesis? I.e. that it estimates the main signal using peak detection and subtracts those from the signal in order to get the noise?

https://en.wikipedia.org/wiki/Spectral_modeling_synthesis

Noise analysis in Spice (Berkeley Spice) is done by summing up the power spectral density from every noise source in the circuit.
There are a couple caveats.

The circuit is assumed to be linear.
In other words the circuit is first solved for a specific DC operating point then each of the components' equivalent linear noise sources are modeled as thermal noise sources.  Each of the components internal resistance values produce thermal noise whose frequency spectrum is shaped by the LINEAR response of the circuit.
Each source is assumed to be uncorrelated with any other.
That means that the power spectral density of the resulting noise is the sum of the power spectral density of each source individually. This is not necessarily the case when you have matched pairs of transistors that are closely coupled and share phonic coupling.
The 1/f noise is a complex parameter model.
1/f noise is difficult to model but there have been several standard models developed. If you stay away from low frequency in the analysis you can avoid dealing with this issue.

See google for more info on noise models such as Noise Sources in Bulk CMOS.


Canonical Rectangular Modulation http://zone.ni.com/images/reference/en-XX/help/371025N-01/loc_eps_rectangular_modulation.gif

I just wanted to make sure that my assumption that in this figure $\phi(t)$ is representing $m(t)$, and $R(t)$ is  the amplitude of the carrier frequency is correct.

Your understanding is not correct. The baseband processing maps the message to a complex-valued signal

$$x(t)=R(t)e^{j\phi(t)}\tag{1}$$

or, equivalently, to an in-phase component $I(t)=\text{Re}\{x(t)\}$, and a quadrature component $Q(t)=\text{Im}\{x(t)\}$. It is not clear from the figure how $R(t)$ and $\phi(t)$ are related to the message signal $m(t)$. The modulated signal is constructed by multiplying $(1)$ by a (complex) carrier, and taking the real part:

$$s(t)=\text{Re}\{x(t)e^{j\omega_ct}\}=R(t)\cos(\omega_ct+\phi(t))\tag{2}$$


Changing the number of cycles of an sinusoidal input signal Imagine we have an ultrasonic sensor that transmits a multiple cycles of sinusoidal signal/waveform that can travel up to a specific distance/range.

Can we achieve more distance coverage using more number of cycles in the input signal (at the same frequency) or this change can not affect the coverage?

I understand that more number of cycles will result in more signal energy. Does that mean that having more energy in the signal will result in less signal deterioration when SNR goes down in longer distances due to attenuation?

So yes, increasing the duration of the transmitted sinusoidal pulse (by increasing the number of cycles transmitted) will provide improved detection at the same SNR.

This is equivalent to a brute force bit duplication used to decrease the bit error rate (hence improve the reliable communication) at a given SNR. Note, however, that you will be transmitting at a slower rate than a shorter window.

So in effect, this method of longer window will compromise transmission rate vs improved detection at longer distance (lower SNR).


How to get around the circular shift property of Discrete Fourier Transform? I understand that when we introduce a linear time shift using DFT on a finite sequence, the algorithm assumes that the signal repeats itself outside of the given range. Here is an example explaining this circular shift property of DFT (Oppenheim, 1998):

Now say if I sampled a finite length (X samples) of a pure sinusoid, but the sinusoid itself is not perfectly repeating when joining the beginning and end together. Hence, when I want to introduce a time shift in the frequency domain, the  DFT algorithm will assume that my signal repeats itself every X number of samples. So I get this:

The output signal (red) is not a smooth sinusoid because my original signal does not repeat itself perfectly.

My question is, is there a way for me to apply the time shift to the pure sinusoid without having this discontinuity in the output? That is, can I have another algorithm that continues my signal smoothly for me after the time shift has been applied?

Thank you all for your time!

Sorry No.

Circular shift property of the DFT (or actually the DFS, @robertbristow-johnson will love this!) is established by law; you cannot get away from it using other clever techniques...

May be you can introduce some redundancies (such as long set of samples but short windows on them, i.e., zero padded signals) you can do some tricks. But at a fundamental level, DFT will always perform a circular shift.


Given a filter (coefficients) how can I make it minimum phase? Consider we have a low pass FIR filter (or other type of filter) that means the filter coefficients are given. Now we want to have the same (or almost) frequency response  but a minimum phase filter. Is there a method to derive the minumum phase filter from the given filter coefficients.

One can find the roots of the given filter, is there a way to bring the roots that are outside the unit circle into the unit circle so that the corresponding roots are minimum phase?

I've seen Julius' MATLAB code and I know what it does.

Essentially, given an LTI filter with impulse response, $h[n]$, and frequency response:

$$\begin{align}
 H(e^{j \omega}) &\triangleq \Big| H(e^{j \omega}) \Big| \, e^{j \phi(\omega) } \\
&= \sum\limits_{n=-\infty}^{\infty} h[n] \, e^{-j \omega n}
\end{align}$$

Then $\Big| H(e^{j \omega}) \Big| > 0$ is the magnitude response and $\phi(\omega)$ is the phase response.

To be a minimum-phase filter, the phase response (expressed in radians) is the negative of the Hilbert Transform of the natural log of the magnitude response.

The natural log of the magnitude response, $\log \big| H(e^{j \omega}) \big|$, is just like decibels (dB) but is expressed in a different dimensionless unit, the neper.  8.685889638 dB is equal to 1 neper.  Essentially the nepers is the real part of the complex natural log of the complex frequency response and the phase in radians is the imaginary part.

Just like radians is the mathematically natural unit for angle, so also are nepers the mathematically natural unit for relative change of magnitude.  A change of magnitude of 1% is nearly the same as 0.01 neper.  But, like dB, going up 0.1 neper followed by going down 0.1 neper will land you at exactly the original magnitude.  But going up 10% followed by going down 10% will land you at slightly less than the original magnitude.

So, the (unwrapped) minimum phase of a filter, given it's magnitude response is:

$$\begin{align}
 \phi(\omega) &= - \mathscr{H} \Big\{ \log \big| H(e^{j \omega}) \big| \Big\} \\
\\
 &= - \int\limits_{-\pi}^{\pi} \log \Big| H(e^{j \theta}) \Big| \, \Big( 2 \pi \tan \big(\tfrac{\omega-\theta}{2} \big) \Big)^{-1} \, \mathrm{d}\theta \\
\end{align}$$

Where $\mathscr{H} \big\{ \cdot \big\}$ is the Hilbert Transform.

If we were to evaluate that integral, we would need to do something called "take the Principal Value" to deal with the division by zero when $\theta = \omega$.  But we won't do the Hilbert transform that way.

Okay, so Step 2 is understanding that the Discrete-Time Fourier Transform (DTFT) of this sequence:

$$ g[n] \triangleq \begin{cases}
0 \qquad & n<0 \\
1 \qquad & n=0 \\
2 \qquad & n>0 \\
\end{cases} $$

is

$$\begin{align}
 G(e^{j\omega}) &= \sum\limits_{n=-\infty}^{\infty} g[n] \, e^{-j \omega n} \\
\\
 &= \frac{e^{j\omega}+1}{e^{j\omega}-1} \\
\\
 &= \frac{e^{j\omega/2}+e^{-j\omega/2}}{e^{j\omega/2}-e^{-j\omega/2}} \\
\\
 &= \frac{2\cos(\omega/2)}{2j\sin(\omega/2)} \\
\\
 &= -j \frac{1}{\tan(\omega/2)} \\
\end{align}$$

i'm running outa time.  i gotta return to this.


Z Transform - Do i always need 2 poles for every "peak"? I am quite new to digital signal processing (and also Z transforms). I am reading about frequency response modelling and have some questions

do we always need a pair of poles for each "peak" in the frequency response ?
why is this so ?
what happens if i dont have a pair of pole for each peak ?

If your system has a real impulse response, then pole & zeros will be in complex-conjugate pairs; so a pole at $\omega = \pi/3$ implies another pole at $\omega = -\pi/3$ .

Further, by the definition of poles, as the frequency $\omega$ get closer to the pole frequency, the frequency response magnitude increases (so as to say peaks!). So for every peak in the (positive frequency axis) then you would have at least two poles at the same frequency (one for positive frequency and one for the negative frequency). Note that there could be more, but at least two is required for a real system.

However note that, having $2N$ poles does not necessarily mean that the frequency response curve has $N$ peaks (sharp rising transitions). Indeed counter-examples are easiy to find, such as an Butterworth or elliptic IIR lowpass filters having no peaks, compared to a sharp isolated rising transition.

Note that in those filters the poles are bundled together such that their indivual peaks have overlapping effects on the total frequency response (yielding a flat passband rather than isolated peaks).


PNG IDAT specification I'm going over the W3 PNG specification (creating a PNG library from scratch) and I've finally found how to create a green 1x1 image.

Now I'm trying to create a bigger image of mixed Red, Green and Blue pixels. Let's say a 4x4 image. Sadly I'm getting all the pixels mixed and some of them are Black or Pink.

Details:

Signature: OK
IHDR: OK

Width: 4
Heigh: 4
Bit depth: 8
Color: 2
Filter: 0
Compression: 0
Interlace: 0

IDATA: 4 chunks*
IEND: OK

*IDATA chunks:

4 bytes length: Number of bytes of the DATA after Zlib deflate
4*1 byte type: IDAT
X bytes data: 4*4bytes unsigned integers, Zlib deflated, one after the other, Network byte order

255 for blue**
65280 for green**
16711680 for red**

4 bytes crc: OK

**Bitwise results of:

alpha<<24 | red<<16 | green<<8 | blue

With alpha, blue, green and red taking values from 0 to 255

What's wrong with this?

I just answered this in Stack Overflow.  Here's the answer again:

The PNG pixels are in RGBA order, not ARGB, so you'd need to write

alpha | red<<24 | green << 16 | blue <<8

But you are writing color type 2, so your pixels should be 3 bytes each instead of four; you can't really encode them in 4-byte integers as you've done. So either change the color type to 6 or switch to writing each sample as an individual byte.

Also you need a filter byte at the beginning of each row. Zero will work for you. So in your example 4x4 RGB image you will need to write 13 bytes per row

0 R G B R G B R G B R G B
0 R G B R G B R G B R G B
0 R G B R G B R G B R G B
0 R G B R G B R G B R G B

and zlib-compress that. It can all go into one IDAT chunk. If you need to write smaller IDAT chunks, you have to zlib-compress the image first, then split the zlib output into pieces that you put in consecutive IDAT chunks.