How to test I/Q modulation with only one branch (I or Q)? Given is a direct-conversion I/Q up- and downconverter system. Receiver and transmitter share the same (10MHz) reference and hence the LO frequency is identical and there is an (unknown) phase difference between them. The transmitted RF signal is given as:

$$
x_{\rm rf} = x_i \cos(\omega_c t) + x_q \sin(\omega_c t)
$$

Now I look at only one channel (e.g. the I channel) at the receiver (that means that any I/Q imbalance of the receiver is irrelevant). At the I branch, this signal is downconverted with an arbitrary phase shift phi and filtered with H:

$$
y_i = H( x_{\rm rf} \cos(\omega_c t + \phi) ) = \cdots \\
\approx H( \cos(\phi) x_i/2 - \sin(\phi) x_q/2 ) \\
= \frac{\cos(\phi)}{2}H(x_i + x_q \tan \phi) \\
= G(x_i + \alpha x_q) \\
= G(x_i) + \alpha G( x_q) \\
= G_1(x_i) + G_2(x_q)
$$

The filter H is contains all filtering effects (mainly anti-aliasing filter).

I want to find a linear relationship between the time-domain samples $y_i[n]$, $x_i[n]$ and $x_q[n]$. A straight forward way is to bring this equation in a vector-/matrix form, model $G_1$ and $G_2$ as FIR filters with $Q$ taps and solve it via Least Squares:

$$
\mathbf{y}_i = \mathbf{X}_i \mathbf{g}_1 + \mathbf{X}_q \mathbf{g}_2
$$

This works nicely and I get the result as expected (in my case, the error between $\mathbf{y}_i$ and the received samples is -45dB). The caveat is I get a set of $2Q$ (real-valued) coefficients but according my derivation, $G_1$ should identical to $G_2$ up to a constant factor. If I would do the same for the $Q$ channel I get again $2Q$ coefficients and a total of $4Q$ real-valued coefficients if I were to model the I/Q relationship that way.

Clearly I can model everything directly in complex domain (with $Q$ complex-valued or $2Q$ real-valued coefficients) but then I cannot look at only the I channel without effect from the Q channel.

Hence a more proper way is to look at this relationship:

$$
\mathbf{y}_i = \mathbf{X}_i \mathbf{g}_1 + \mathbf{X}_q \alpha \mathbf{g}_1
$$

However, this is not a least squares problem any more. For that reason, I sweep $\alpha$ within a reasonable range, compute the LS solution and take the $\alpha$ with the minimum error (this seems to be convex). But if I do this, the error -35 dB (as opposed to -45dB above).

I really do not understand why this approach would not work.
Why? And what is the proper way to model this?

Background info: I am testing my own direct conversion receiver. If I transmit data in only one channel (I or Q) and look at only one output channel (e.g. I) I get the performance that I expect based on linearity-, noise- and phase noise measurements of the receiver. If I look at I/Q at the same time (and model the the effect of the receiver as complex-valued FIR filter) then the error is more than 5 dB worse. For debugging I would like to avoid any interaction (such as I/Q imbalance) between I and Q and only look at the I channel output. But as soon as I transmit a full complex-valued input and look only at the I output my result is again at least 5 dB worse than predicted by noise, linearity, phase noise.

As another step, I use an Arbitrary Waveform Generator (Rohde & Schwarz SMW200A) used as transmitter and a VSA (Rohde & Schwarz FSW) as a receiver where I see the same issues. Both instruments are calibrated and the transmitter has superior image rejection. For these reasons I assume I made a conceptional or algebraic mistake in my derivations/assumptions.

The FSW is superheterodyne.

This means that it will not have a DC offset or I/Q imbalance, but the IF path will appear as an asymmetrical lowpass (because it is really a shifted bandpass), which is then compensated for in the calibration.

The first question: have you run the calibration properly?

leave both machines running for half an hour in a stable environment.
disable the generator output (OUTP OFF)
calibrate both (*CAL)

Lots of equipment lenders don't understand that the calibration sequence is used to compensate temperature, run the calibration once (usually on a cold machine) to make the error message go away and tell customers to never touch that button. This gives a fairly bad filter to compensate the IF, and since that filter is applied before downconverting the IF to zero, this will also appear as a complex FIR filter to you.

Second, test your setup with a simple test:

generate a sine wave (*RST, FREQ 2 GHZ, OUTP ON)
look at the vector display in IQ Analyzer (*RST, INST IQ, LAY:REPL '1',VECT)

You should see a single stationary dot there, with a position that corresponds to an arbitrary phase angle that is the sum of the phase between the LOs (they are certainly not the same frequency!) and the digital downconversion inside the FSW, which together give a constant, but only as long as all the oscillators keep running at their fixed frequency.

Once you have reached that state, you should be able to repeat your measurement with more complex signals.


What is the expression for intensity of Laser beam focussed at a point at an angle theta? If the laser beam is focussed directly at a point then the simplified expression is:
$$I_1(i,j) = 8\frac{P}{\pi d^2} e^{-\frac{8}{d^2}\left(x^2(i)+y^2(j)\right)}$$

If the laser is focused at same point but directed at an angle theta then what will be the intensity expression?

I am asked to work on a project related to quadrant detector. In this I have to plot the current variations of the four photodiodes. For finding the current I need the intensity of laser beam reaching the four photodiodes. So in case I focus the laser beam directly at a point then intensity is same in all four quadrants and that is given by expression above. If I focus the laser beam at an angle theta then what will be the expression of intensity. I want to know the intensity expression, so that I can plot it in Matlab.

I tried to replace d in above expression with d/cos(theta) but that is not correct. Next, I replaced x(i)^2+y(j)^2 with (x(i)^2+y(j)^2)*cos(theta)^2 but that is also not correct.

Please help me out with the correct result.

Methodology
We can apply simple ray optics to this problem, and it gets easy!
What you're describing is an Orthographic Projection (OP).
So, what you'll do is

realize circles subject to OP become ellipses
take the center of the circles of equal intensity (in your case, $x^2+y^2=\text{const.}$), and project that. It'll be the center of all the ellipses of equal intensity.
figure out what the major and minor axes of the ellipses are going to be
simply write down the formula for points on ellipses with given axes.

Orthogonal Projection
For the math behind OP itself, I'd like to refer you to the wikipedia link given above. It's really just a single matrix multiplication, if your original (circular) image plane and the projected plane both contain the origin $(x,y,z)=(0,0,0)$; it becomes especially easy if the intersection of these two planes is e.g. the $y$-coordinate axis, which we can assume here; you can calculate the normal vector of the slanted plane simply by multiplying the cosine (or sine, depending on your construction) of the incident angle to unit vector on the $z$ axis, then.
Major axis of the ellipses of equal intensity
For the question of what the major axis is going to be: it's the projection of the diameter of the original circle that's parallel to the plane intersection (source). Luckily for you, you've defined that intersection to be a diameter of the circle! Since all points on the intersection of original and image plane are invariant to the OP, you have to do nothing to find the major axis of any circle of given intensity: it's simply the distance on the $y$-axis of the original!
Minor axis
The minor axis is simply perpendicular. You'll find that the projection is also especially easy, and is simply the "complementary" trigonometric function to what you've used to find the normal of the projection plane.
Overall formula
Now that we know how to convert perpendicular incidence to slanted incidence circles, we just have to figure out how to do the inverse. Luckily, since the projection is such a well-formed matrix, that's just taking the inverse– and that's trivial, given the structure of the matrix!


Quadrature demodulation, homodyne detection, lock-in detection - what's the difference? Quadrature demodulation

Homodyne detection

Lock-in detection

What is the difference between these terms? Are they all synonomous? In my mind they kind of are, so I'm never sure which to use.

I became familiar with this concept first through the use of lab devices called lock in amplifiers. However, homodyne detection seems to be the same thing. Quadradture amplitude demodulation seems to be the best term, since it describes what you're actually doing (I have no idea what the word 'homodyne' even means).

They are essentially the same but may have key differences varying across implementation. The process is identical in terms of multiplying a signal by a sine wave at the center frequency of the signal (the carrier) which moves it to DC, and then low pass filtering that signal to maximize the SNR for the detection (presence of the signal at the carrier).

For the differences that may exist the homodyne and lock-in amplifier may or may not be done with quadrature demodulation. Note the following from the relationship of multiplying two sinusoids:

$$\cos(\alpha)\cos(\beta) = \tfrac{1}{2}\cos(\alpha+\beta) + \tfrac{1}{2}\cos(\alpha-\beta)$$

If we multiply two signals of the same frequency in phase and low pass filter to remove the higher frequency term (as denoted by the LPF[]) we get:

$$\mathrm{LPF}\big\{\cos(\omega t)\cos(\omega t + 0°)\big\} = \tfrac{1}{2}\cos(0) =\tfrac{1}{2}$$

but if in quadrature we would get:

$$\mathrm{LPF}\big\{\cos(\omega t)\cos(\omega t + 90°)\big\} = \tfrac{1}{2}\cos(90°) = 0$$

So a key question is what are you trying to detect? When trying to detect phase, the quadrature demodulator would actually be best as it would have the highest sensitivity to phase changes (the slope of the output voltage vs phase shift of the two signals being multiplied is highest when the signals are in quadrature). Similarly it will detect with maximum amplitude the presence of a signal that is also in quadrature (as this would be the first case essentially if both signals are at 90°). For arbitrary detection of any angle you would use an IQ demodulator that would multiply the input signal by cosine and sine with I and Q outputs each low pass filtered to show the detection for any phase of the input.


How to establish transfer function of a speaker? I'm working on a signal processing project about vocal system, and I'm trying to use controlling theories to solve the problem.

I need to get the transfer function, $H(s)$, of a speaker, from electrical signal to voice. But it seems that I have to use a microphone to pick up the voice. Considering the microphone itself has a transfer function $G(s)$, I will finally get $G(s)H(s)$ rather than $H(s)$.

So, is there any method to get TF $H(s)$ from electrical signal to vocal signal?

Well I have another problem. I've seen someone else estimating TF $G(s)H(s)$ by experiment. The bode diagram looks so strange and complex that I think it will be difficult to get a mathematical description of the system. So what method is usually used in controller-design when the object has a complex freq. response?

You can use a spark gap to measure the impulse response of the microphone. Spark gaps discharge within a very short time interval, producing a very broad spectrum that is effectively constant in the audible range.

If you also make sure that reflections that reach the microphone are arriving much later compared to the duration of the impulse response, you can separate the room response very well from the microphone response. Repeated measurements improve the estimation.

With the microphone response known, you can solve for the speaker response. However, technically you still have the room response and directional dependence of the speaker and microphone that you might want to consider.


Identifying the boundary points of a shape based on the centroid I have a set of binary 2D images with objects in white and their boundaries in black (attached example). I also have a list of centroids (x,y coordinate pairs) identifying all objects.

For each object, I am trying to identify all the x,y coordinates of that object's boundary (black pixels). Anyone have an idea on programming logic that can help achieve this?

I apologize in advance if this post belongs in another section. If anything is unclear I can clarify. Thanks.

In general you would be looking at "Image Segmentation" techniques.

Since you do have the centroids, you can use them in "region growing".

In this technique, you start from a seed point and depending on the sort of pixel connectivity you are interested in you expand the identified region. For example, with 4 connectivity, you would start at the seed point, which is clearly part of the region. You would then examine the 4 immediate neighbours of the seed, each one having 2-3 neighbours (there will be overlaps). If those pixels belong to the region (e.g. they turn out to be white), you examine their neighbours and the neighbours neigbhours and so on.

The algorithm terminates when there are no more neighbours to consider for their connectivity.

You can write this from first principles, or you can use skimage's label(). In fact, that function applies pixel connectivity slightly differently, by building regions up as it reads the image, creating a new one if it finds a disconnected pixel, so, you would not even need the centroids in that case.

Hope this helps.


Electric Bass with RJ45 Digital Output I'm working on a project bass guitar with individual outputs per string using Nu pickups from Cycfi Research: https://www.cycfi.com/
The common sense approach is to simply take their circular connector and use this to interface, but I'd like to know how realistically could I use a multichannel ADC/MCU to convert the 4 (or more with scalability accounting for 5 strings or two pickups) outputs serially to an RJ45 in real-time?  The goal would be to preserve each individual output to process digitally in a box or use a DAC to split the channels back out to process with analog pedals (for example).
This is a question of "could" not "should," and it comes from a place of learning, tinkering, desire to use a more commonly available cable for lower cost (the circular is a good deal for its quality, but it's $84 and not common), and the ability to send power all in one cable.  Development time, headache, and cost of electronics inside or outside the guitar are not a factor, this is a hypothetical since I just don't have experience working with multichannel ADCs or ethernet and I'd like to use this as an alternative project to learn from.
I don't mind if this is wrong, I'd just like to understand if this is a rabbit hole or the elusive unobtainium.  Thanks in advance!

I'd just like to understand if this is a rabbit hole or the elusive unobtainium.

A little bit of both. Ethernet is really NOT a great choice for real time audio. You would need a fairly beefy processor that can handle a network stack (TCP or UDP). Getting the latency down to something that's playable is extremely challenging.
Even if you get it on Ethernet, where would it go ? There some multi-channel network formats but these tend to be in the PRO space  and expensive.
If you just care about the connector, you can certainly run I2S or TDM through an RJ45 + Ethernet cable. That works quite well (and I've actually managed to incorporate this into a commercial product). Again, you would need some type of receiver that's matched to your transmitter.


Does low-pass filtering a square wave cause a 'ringing' to develop in the signal? Using the figure below as an example, I wonder if it is possible that the small oscillations in the signal at the time of each step shift could be caused solely by the application of a low-pass filter to a pure square wave?

These ringing artifacts are due to the Gibbs phenomenon. They are indeed caused by lowpass filtering a function with discontinuities. What you see is basically the oscillations of the step response of the lowpass filter. The oscillations can be reduced (or completely eliminated) by choosing a lowpass filter with less or no ringing in its step response. E.g. a Gaussian filter has no ringing at all. However, lowpass filters with well-behaved step responses perform worse in the frequency domain in terms of passband deviation, transition band roll-off, and stopband attenuation.


Realize Simulink block "Ideal ADC Quantizer" in Matlab Yesterday I posted following topic: "FFT alternates '0' for a quantized SIN wave".
During post editing I was perturbed by some affair, I didn't formulate correctly the problem. So, my post was blocked (i.e. "hold on"). My subsequent update didn't modify the post status: it remains "hold on".

Since, I advanced in the problem resolution. The problem originates from erroneous quantization of SIN waveform. For quantization I use Matlab "quantiz" function, that take 2 parameters: "partition" and "codebook". After changing the fragment of code that is related to quantization (please, see below) the problem disappears, i.e. FFT displays correctly for all frequencies of SIN waveform.

dy = (2*A)/2^res;
partition = (-A+3/2^(res+1):dy:A-3/2^(res+1));
codebook = (0:2^res-1);
[index,yq] = quantiz(y,partition,codebook);
yq = 2*yq/(2^res-1)-1;

where:

A - amplitude of SIN wave (=1)

res - resolution of quantizer (10bit)

But nevertheless some incertitude remains: when I to solve the same task in Simulink (using "Ideal ADC Quantizer" block, the results of Simulink simulation are slightly different from Matlab simulation.

Any suggestions ?

Thanks

P.S. Here is Simulink model (with ADC settings)

Here is also comparison of 2 FFTs:

SIN waveform, quantized in Simulink (TOP)
SIN waveform, quantized in Matlab and imported in Simulink (Bottom)

Resolved.
The problem was in erroneous quantization.
Here is the correct quantization code:

% Quantization
dy = (2*A)/(2^res-1);
partition = (-A+A/2^res:dy:A-A/2^res);
codebook = (0:2^res-1);
[index,yq] = quantiz(y,partition,codebook);
yq = 2*A*yq/(2^res-1)-A;

where:

y - SIN waveform values

A - amplitude of SIN waveform

res - resolution of quantizer


Garment Cropping from mannequin I have two images – mannequin with and without garment.

Please refer sample images below. Ignore the jewels, footwear on the mannequin, imagine the second mannequin has only dress.

I want to extract only the garment from the two images for further processing.

The complexity is that there is slight displacement in the position of camera when taking the two pictures. Due to this simple subtraction to generate the garment mask will not work.

Can anyone tell me how to handle it?

I think I need to do registration between the two images so that I can extract only the garment from the image?

Any references to blogs, articles and codes is highly appreciated.

--
Thanks

The approach you are alluding to involves a) first registering the images and b) using a subtraction operation to remove the background (which is the mannequin).

First on (a) (registration): Image registration deals with figuring out the displacement of corresponding pixels in two images. The problem with real-world images is that the world is 3D and because of the depth in the scene, nearer surfaces will move more than surfaces deeper in the scene. Also, a surface can self-occlude certain parts of the scene. I only state these so that you are aware of the complexities involved. Having said that, in the scenario you have posted above, you might be able to do very well with a simple translation model for the displacement. In other words, you can assume that all pixel in the second image move by $(\delta{x}, \delta{y})$ from its original position. So if a pixel was located at $(x,y)$ in $I_1$, then its location in $I_2$ will be $(x+\delta{x}, y+\delta{y})$. This can be assumed to be true for every pixel (global translation). You now need to estimate the horizontal shift $\delta_x$ and the vertical shift $\delta_y$. There are many approaches to this problem. Here is an academic reference to registration methods (may be TMI) Zitova and Flusser. You may wish to start off by trying the ITK toolkit and try using the "translation transform". If you find that after accounting for a translation transform, the image subtraction still leaves parts of the mannequins, you might want to try some more advanced transforms. For example, the camera may have rotated slightly in your hand. So in addition to a shift, you can also introduce a rotation. Also, you might have taken one shot closer to the other, so there may be a difference in scale. This necessitates the use of an affine transform which is also included in ITK. Barring these "global" methods which assumes that every pixel moved in the same way, you can also use local methods. For example, you can estimate a $(\delta_x, \delta_y)$ for every patch in the image. You can divide up the image in patches of $7 \times 7$ pixels and search for the best corresponding patch in the other image. Every pixel within this $7 \times 7$ window is assumed to have moved in the same manner as the motion that was estimated for the patch. A last resort, but the most general is to compute dense optical flow for every pixel in the image. This gives you a shift for every pixel in the image. You can find an optical flow implementation in the popular OpenCV software.

Finally, on (b) (occlusion): You will inevitably end up with some pixels without a valid shift because they simply aren't present on the other image. Or some pixels are uncovered in the second image which were covered on the first. See these paper for clues Occlusion aware optical flow and TV-L1 Optical Flow with Occlusion Detection


Sequentially get the outer contour of an Edge Image Hi CV/Pattern Recognition Community,

because my last thread was concerning the segmentation as BLOB and started to getting too big, with no clear results.
I would like to reconcider the last step of that process regarding the sequential gaining of object contour points.

I need to extract the contour of an image from a defined starting point (left) to the right (endpoint).
As you can see in the next image, there are two yellow windows (left and right). These are the mentioned start- and endpoints.
And what I want to do: Starting from the left yellow window and moving to the right, saving (x,y) coords of the regarding points. This is what the red line (contour) does.
The sequence of the contour pixels does matter.

I am at the step where I've only got a binary image,
So far, so good.
In the worst case still with holes and possible noise in the rest of the image.

So, here are some images, which I want to extract.

For this example, I'd like to get the outer contour line of the object on the table/underground. As you can see in the right upper corner and in the left upper corner are fragments of the edge detection process, which should not be traversed.

Thank you in advance.
Matlab and/or OpenCV Hints are welcome!

You can do the following in order to get a connected contour. Get a few points on the contour say like 100 points or so. And then fit a cubic-spline on it.

I had done something which might be helpful to you. Below is what I had done.
The project was clustering of cars by shapes. So I had images of cars and I a user had marked the boundaries of the car. Red dots are the manually marked points. Then I had fit a closed cubic splines on these few points. The blue contour is the cubic spline. I then used this contour for my purpose

I am sure you can do something similar to get the contours. Please be aware that it is possible to get an open curve with cubic splines. Let me also tell you constructing splines from these set of co-ordinate points is simple. What you basically want to do is fit a cubic curve between each set of points.

Please ask if you need further clarification on what I think.


where is the matched filter being used exactly in Rx PHY? Recently I learned about the matched filter and its usefulness as it maximizes SNR at time t0.
So I looked up online to see how this amazing filter is being used exactly in 802.11a Rx PHY, but couldn't really point out which blocks in Rx PHY the matched filters are used exactly.

I suspect that matched filter can be used as a demodulator, but if one deploys 64-QAM or 128-QAM, do we still use matched filter for the demodulator in this case as wel??

so I have two questions :
Q1: Where in Rx PHY is the matched filter being used exactly?
Q2: If matched filters are used in demodulator as I expect, how is it implemented for a complex demodulation scheme (i.e., 64-QAM).

Q1: Where in Rx PHY is the matched filter being used exactly

Not at all / it is a rectangular filter.

802.11a is an OFDM system, and that means it's a multi-channel thing – it divides the wide channel into many narrow channels, each of which is flat, i.e. only a single complex coefficient. It uses a rectangular pulse shape on each of these.

Thus, there's no matched filter explicit anywhere in the receiver.

When you further your study of modulation techniques, you'll learn why OFDM has to use a (in time domain) rectangular filter for each subcarrier (because that achieves the zeros in frequency domain that make the "O" in OFDM).

  Q2: If matched filters are used in demodulator as I expect, how is it implemented for a complex demodulation scheme (i.e., 64-QAM).

Filtering is not part of the symbol decider, but comes before it (in systems that, unlike OFDM, do make use of it), after the downmixing (if applicable).


Phase of a real number Could someone please explain in what case the phase of a real number is equal to -pi (and not pi)?

I know that for positive numbers, the phase is zero. For zero, we define the phase as zero as well. And for negative numbers, the phase would be pi. But I was reading some script and there it says the phase of a real number is either 0, pi, or -pi.

Or $2\pi$, or $3\pi$, or any integer multiple of $\pi$.  Any odd multiple corresponds to -1 + 0i and any even multiple corresponds to 1 + 0i, aka -1 and 1.

"Phase of a real number" is a little bit of a misleading label.  What is required here is an understanding of the complex plane and what "phase" means in terms of a DFT bin value.

Your question is equivalent to "For what values of arg(z) is z a real number?"

If that is meaningless to you, I suggest you start by reading two blog articles of mine:

The Exponential Nature of the Complex Unit Circle

And the newest:

Angle Addition Formulas from Euler's Formula

There are of course many other searches.  Your terms should be "complex plane real values" for a start.

This is essential foundation material for a lot of DSP concepts.


Carrier and timing not properly locking? I have designed a simple BPSK system which is operating in UHF band. My transmit data rate is chosen to be 500 kbps. My receiver has AGC, timing recovery and carrier recovery bocks.

I am using Gardner timing recovery with 2 samples/symbol. The loop bandwidth for 2nd order loop filter for timing and carrier recovery is chosen to be

My signal level varies from -30 to -10 dBm at the receiver. I am observing the lock in carrier and timing is not constant. Is it because of the loop parameters or signal level or the header bits chosen for the carrier and timing recovery?

I recommend that the OP confirm via testing each of the loop parameters separately in terms of convergence time and lock-in range. Step response testing is a good approach for convergence time and the range of the step can be increased to both evaluate lock-in range and close-in "linear" response that would affect tracking versus further-out "non-linear" response that would affect acquisition.

To do these one at a time I suggest for example the AGC: with the system locked in carrier and timing, introduce a very small step in amplitude and evaluate the time domain response and how long it takes to converge. Increase the range of this step and the response time should stay the same when in linear operation, and then start to increase or change as the operation becomes non-linear (if it does), and then eventually lose lock if the step is outside of the acquisition range (or more likely causes the carrier and timing to lose lock).  Ideally this is done in simulation such that the other parameters can be forced to remain locked.

Carrier and Timing would follow the same approach: for carrier the level should remain constant and carrier and timing acquired and then introduce a small carrier offset, etc...


Training a NN classifier on a single\double channel out of a surround sound dataset I want to train a neural network for a classification task on the VSD2014 dataset. I have downloaded the movies, but they have a 6 channel audio format (surround sound). 6 channels will cost a very long training time, and so, I wish to trail on a single\double channel first, in order to prove the concept.
Is it better to convert into a stereo or single channel? If so, how do I convert this on MATLAB? Can I get away with, simply, choosing one of the channels? If so, which one? I am guessing the subwoofer channel is not the best choice. I do not even know which one is the subwoofer channel.
How do I train a NN classifier on a single\double channel out of a surround sound dataset?

Generally the order of a 5.1 channel signal is as follows: L (Left), R (Right), C (Center), LFE (Low frequency enhancement), Ls (Left surround), Rs (Right surround).
What you need is a downmix and there are many downmixing methods. Here's the recommendation of ITU:
$$
L' = L+0.7071C+0.7071Ls
$$$$
R' = R+0.7071C+0.7071Rs
$$
However, multichannel signal may contain some hidden information for your classification task. To get a better performance, I don't think downmixing to single channel is a good idea. Stereo signals may be a compromise between performance and training time.


Why does the specgram of this simple sine wave look strange? I'm trying to test some intuition using the matplotlib specgram routine. I first generate a simple sine wave

import matplotlib.mlab as mlab
import matplotlib.pyplot as plt
import numpy as np

t = linspace(0,100,1000)
omega = 2*np.pi*1
y = np.sin(omega*t)
plt.plot(t,y)
plt.show()

This generates the expected image of a sine wave

I can also look at the mod squared of the FFT of this signal:

the_fft = np.fft(y)
spacing = t[1] - t[0]
ps = the_fft * np.conj(the_fft)
freqs = np.fft.fftfreq(len(t),spacing)
idx = np.argsort(freqs)
plt.plot(freqs[idx],ps[idx])
plt.show()

which gives me the square of the Fourier transform, which has a peak at 1 as expected (since my frequency is 1).

I then try to compute the specgram, as

_ = plt.specgram(the_fft,NFFT=100,noverlap=0,Fs=1./spacing,window=mlab.window_none)
plt.show()

which then gives me this .

From this it's not easy to tell that the only frequency present in the signal is 1. What needs to be done in the specgram call to see this? Of course this is probably the expected output of this, so I'm really also asking a basic DSP question at the same time.

You are computing the spectrogram of the DFT as if it were a time signal. Try this instead

import matplotlib.mlab as mlab
import matplotlib.pyplot as plt
import numpy as np

t = np.linspace(0,100,1000)
omega = 2*np.pi*1
y = np.sin(omega*t)

spacing = t[1] - t[0]
plt.specgram(y,NFFT=100,noverlap=0,Fs=1./spacing,window=mlab.window_none)
plt.show()


Sampling rate and number of samples I am doing an experiment where i have to excite a plate and scan a region on it using a laser signal.
I am acquiring the signal from the laser continuously using an NI daq tool. i then have to divide this  signal into a set of discrete intervals, (each corresponding to a 1 mm segment in the scanning region), in which the time between the segments should have integer number of harmonic excitation wave cycles.
The speed that the laser moves in on the scanning region is 1 mm/min.
The excitation frequency is 100Khz, and the sampling rate of the daq tool is 1Msps.
however, i am not sure what should be the number of samples N for the tool?

Sampling Rate of the daq tool tells you how many measurement samples will be collected per second. So, if the sampling rate $f_s=1Msps$, it means that your daq tool is collecting $10^6$ samples per second.


How we do frequency shaping and amplitude shaping in hearing aid My project name is Digital Hearing Aid with frequency shaping and amplitude shaping using adaptive noise cancellation DSKC6713 DSP kit .
Almost I have done this noise cancellation using code studio ver 3.1( specifically used for this kit)
Now I get stuck on frequency shaping and amplitude shaping ....after so many searches I still dont get any coding or algo ...can anyone help me..??

Merium:

To develop a better understanding of how frequency shaping is used in hearing aid processing, consider the following paper, "Physiological assessment of contrast-enhancing frequency shaping and multiband compression in hearing aids", Ian C. Bruce, Physiological Measurement, Volume 25, pages 945-956, 2004.  It is an older paper, but it provides a good overview of some ideas you may want to consider using in your project.  If you are a registered student at a college or university, you should have no trouble finding this article using the library's resources.

Due to the fact that common hearing impairment is associated with loss of high frequency sensitivity, I would suggest implementing a simple filter that provides gain for high frequencies and maintains the gain at lower frequencies.  You could do this using a multiband approach, which would give you more flexibility in setting the gain for each band, or a across single band that might be more readily realized in the DSKC6713 framework.  Have a look at this article to understand the different types of hearing loss.

As far as implementation goes, I would suggest using Matlab to develop and test your processing scheme prior to implementing it in the DSP C++ framework.

I hope this helps.

Michael.


Compute dynamic range in a linear quantization system I know the dynamic range is the ratio between the maximum amplitude and the minimum amplitude of a signal, calculating the dynamic range in dB I get a difference of logarithms:

But why the dynamic range of a N-bit digital system with linear quantization is 6N dB?

Just to summarize the comments, for the next person who comes looking:

decibels are $20 \log_{10}{x}$ (or $10 \log_{10}{\left( x^2 \right) }$ when power levels used instead of the amplitude levels)

When looking at bit noise limited N-bit system, there are $2^{N}$ measurement levels, and the smallest distinguishable difference is $1$. So

DR$_{dB}$ = $20 \log_{10}{\frac{2^N}{1} } = N*20\log_{10}{2}\approx 6 N$


Should I pass Kalman Filter absolute or offset-from-mean sensor values? I'm using Kalman filters to segment the loudness of an acoustic signal from surrounding noise. The problem I've encountered is that muffled or faulty microphones measuring 'silence' (-70dB, -69dB, -71dB, ...) have low variance and so the Kalman filter relies on them as being more stable, "pulling down" the true value.

To solve this, I had the idea to pass in the offset-from-mean of all the signals, which would increase the variance of the "dead sensors" and I could add back the mean to the state estimate delta. But then it occurred to me that maybe the Kalman filter is already doing this by minimising "least squares error"?

A secondary problem I have is dynamically adding or removing sensors without losing my covariance matrix. Currently when sensors leave, I keep the same matrix, but I set it's observation variance R_ij = Infinity, which seems to work. However, I would love to know how to re-arrange my matrix values when adding or removing additional sensors.

As mentioned in another answer, the variance of the measurement noise is a property of the model which is baked in to the Kalman filter.
I think you're running into a different problem, which is that your model is different from reality.  It sounds like you're modeling the system assuming that all the microphones are 'live'.  If that's the case, then the filter won't take a dead microphone as something to be ignored -- it'll take a dead microphone as a live microphone that's 'hearing' utterly valid silence.
Note that Kalman filters can do some apparently absurd things in the presence of a model mismatch.  Because they assume Gaussian distributions for everything, and because the PDF of a Gaussian drops off so profoundly as you stray from the mean, they really don't do a good job handling real-world processes that have low kurtosis, nor do they handle model mismatch very well.
You need to modify your filter, or sanitize its input.  It's probably best to modify the filter by augmenting it with an external "dead microphone" detector (and make sure you get it right -- you're basically automatically throwing out outliers, and there's known problems with doing that too enthusiastically).  Or, make sure a human goes through the data and makes sure the microphones are live.


DFT of complex signal using Goertzel algorithm in C I am implementing a BFSK frequency hopping system with TX and Rx modules. I am using Goertzel Algorithm at the receiver end to demodulate the data i.e. to determine the carrier frequency of the recevied signal data. Following is the implementation of the goertzel implementation:

float goertzel(int numSamples,int TARGET_FREQUENCY,int SAMPLING_RATE, float* modData)
{
int k,i;
float   floatnumSamples;
float   omega,sine,cosine,coeff,q0,q1,q2,result;

floatnumSamples = (float) numSamples;
k = (int) (0.5 + ((floatnumSamples * TARGET_FREQUENCY) / SAMPLING_RATE));
omega = (2.0 * 3.1416 * k) / floatnumSamples;
sine = sin(omega);
cosine = cos(omega);
coeff = 2.0 * cosine;

q0=0;
q1=0;
q2=0;
result=0;

for(i=0; i<numSamples; i++)
{
    q0 = coeff * q1 - q2 + modData[i];
    q2 = q1;
    q1 = q0;
}

result = sqrtf(q1 * q1 + q2 * q2 - q1 * q2 * coeff);
return result;
}

The implementation works well for a signal where the received signal is a real signal ( a set of float values). Would this algorithm still work for complex signals. The complex signal would be typically be represented by a set of float values, with each alternate element representing a real and imaginary value. In this case, how should the above goertzel function be used, should it be used to calculate the output for the real and imaginary components separately, or should it read the entire complex signal array elements in a sequence (both real and imaginary array elements) and then  calculate the magnitude at the end of it. If not, can some please point me to how to get the goertzel output for a complex signal samples. I tried googling it but could not find anything concrete about how to use goertzel algorithm on complex signals.

Also, I wanted to ask one more thing. If my signal has carrier frequencies at both 4000 and -4000 Hz and a sampling frequency of 18000 Hz, then for a real signal as input, the goertzel output is the same for input frequency of 4000 Hz and 14000Hz (this is what should be). But since 4000 Hz and -4000Hz are one of my carrier frequencies, I am not able to determine whether the received signal is at 4000 Hz or -4000 Hz. Is there any method to do this without doing a frequency shift. I am able to differentiate using frequency shifting (multiplying the received signal with a sine wave before demodulating). But I am looking for a method to do this without using the frequency shift technique.

--Anshu

The Goertzel algorithm (which is really just an efficient way of calculating what amounts to a single DFT bin at an arbitrary location) is defined for complex input, just like the DFT. A real input signal is really just a special case where the imaginary part is equal to zero. As far as your last question is concerned, if your input signal is real, then its spectrum is conjugate symmetric. That means that the spectrum value at -4000 Hz is the conjugate of the value at +4000 Hz. If you want to differentiate between the two frequencies, you'll need a complex input signal.


Bilinear transform with pre-warping for systems other than classical filters I only seem to be able to find online information about applying bilinear transform + pre-warping to filters (like butterworth, etc.) with only one edge frequency that is purposely 'designed' into it. So the BT is easy to apply, just substitute with the usual tangent formula including edge frequency and sampling frequency.

But what if my analog rational transfer function describes a physical system (like an elastic mechanical body with several eigen-frequencies, or an electric circuit network)? Do I have to apply the BT to every pole/zero separately then? And what if the transfer function comes from measurements - I don't know the poles/zeros, because I did not 'design' the 'filter'? Does it really mean I have to factorize it (which we know is sensitive to error) just to transform it from analog to digital?

PS, Sidenote: I am going to implement the numerics involved (in C#) myself. So Matlab is certainly a nice thing, but it does not help me if not all relevant code is publicly available.

The bilinear transform ("Tustin's method") is indeed mainly used for transforming frequency selective filters with magnitude responses that are optimal with respect to some criterium, such as Butterworth, Chebyshev, or Cauer filters.

For more general systems, the bilinear transform is usually not the best choice because of the frequency warping, which mainly affects frequency domain behavior close to Nyquist. For frequency-selective filters, cut-off frequencies can be pre-warped such that the cut-off frequencies of the resulting discrete-time system are correct, but for more general systems it is usually not clear which frequencies should be matched and which shouldn't. Also, for filters with approximately piece-wise constant magnitude responses, frequency warping is not a big problem as long as the cut-off frequencies are mapped correctly. For general systems one would usually prefer the error between the continuous-time system and the discrete-time system to be distributed more evenly in the frequency domain, which is not possible with the bilinear transform.

There are several methods for converting continuous-time systems to discrete-time systems. All of them have advantages and disadvantages, and it depends on the application which one is to be preferred. Take a look at Matlab's c2d function, which implements several methods. Also take a look at this question and its answers.

For approximating a measured response, you would directly design the corresponding system in the discrete-time domain. For general frequency responses, the frequency sampling method is very simple and effective. If an IIR approximation is desired, the Matlab function invfreqz could be used. It tries to fit a rational transfer function to given frequency response data.


psd of an ac voltage signal In an lab class experiment who's subject of interest was Johnson's noise I measured the ac voltage of a battery. Now I wish to find the psd of the measurement and show that it is constant.

1.Is such a thing even possible with 255 ac measurements in a $0.043s$ spacing?

How can it be done without a direct measurement of the voltage?

To compute the PSD, you'd need full amplitude/phase measurements of the noise signal.

If by "255 ac measurements", you mean measurements of the RMS amplitude by a digital multimeter, then you do not have the necessary phase information.

If by "255 ac measurements", you mean you are using an oscilloscope that is AC coupled, then you can consider the measurements to be a measurement of the voltage with the average value ("DC component") removed.   If we consider the voltage signal to be
$$
v_\text{total}(t) = v_\text{DC} + v_\text{AC}(t),
$$
$v_\text{DC}$ is going to be nearly entirely the "normal" contribution of the battery, and constant with time, and $v_\text{AC}(t)$ is going to vary with time and due entirely to thermal noise.  If your measurement is measuring $v_\text{AC}(t)$, then that is what you want anyways.

Once you have waveform samples, the PSD can be estimated using the periodogram (easiest, but noisy) or a more advanced method like Welch's or Blackman-Tukey method.


Time delay estimation in low SNR So I am interested in finding delay between two noisy narrowband signals,for example, two arterial pulse waveform. I can easily get a correct estimate of delay using cross-correlation method, but it fails in low SNR cases. I know Generalized Cross Correlation (GCC) method is supposed to give better results in these cases, but as I am working with narrowband signals, and also due to incorrect estimates of power spectral densities, GCC is very unstable, and the results are worse than normal cross-correlation. So my question is are there any methods for finding time delay between signals in low SNR and are more accurate than the above 2 mentioned methods?

Quite often, performing a time-delay estimation in a transformed domain can help,as long as the signal gets better concentrated, and/or the noise is more spread.

So performing the CC or GCC in a time-frequency or time-scale (wavelet) domain could be useful, as long as you can select the frequency  band or the scale in which the correlation will be meaningful.

If the noise characterics are unknown, you can benefit from wavelet whitening capabilities.

You can check for instance:

On Wavelet Denoising and its Applications to Time Delay Estimation, 1999

For Matlab codes:

Wavelet Cross-Correlation for Lead-Lag Analysis
modwtxcorr: Wavelet cross-correlation sequence estimates using the maximal overlap discrete wavelet transform (MODWT)


FFT looking messy after filtering process I have a vibration signal (25.6 kHz)  which has been through the following filters:

1: highpass (2000 Hz) -> (I'm not interesten in the harmonics of the system. I then perform spectral kurtosis to find the optimal bandpass bandwidth.)
 2: Bandpass (8750 Hz,9250 Hz) ->
 3: Rectification (abs(xi)^2) ->
 4.Lowpass (2000 Hz)

The fft of the final signal looks like this:

The Bandpass filter bandwidth and center frequency was selected based on the kurtogram of signal after initial lowpass.

I am no expert in signal processing- does the FFT look "correct"? It seems odd that the amplitudes rarely go below 0.05...

Updates

The original signal (25.6 kHz):

Zoomed in in original signal (25.6 kHz):

Final signal in time domain:

FFT of original signal (before step 1)

Objective of project:

To look for any fault development in the early stages of a gearbox fault on a wind turbine. I have signal data 412 days, in the form of 412 10-second intervals.
Want to extract features from time and freq. domain and cluster the 412 intervals, to see if there is a pattern.

If your description is correct, you are seeing noise.

After step 1, everything above 2 kHz should be gone. After your band-bass everything below 8750 Hz is gone, so you basically end up with a null set. Since your filters are not infinitely steep there is still something non-zero left over but it's mostly going to be noise and very poorly defined.

Rectification is also a highly non-linear operation, so chances are you also going to get some aliasing in there as well.


How would I calculate the state-variables of digital biquad filter (direct form II transposed) at steady-state? I'm attempting to calculate the steady-state state variables for a digital biquad filter direct form II (transposed). Illustration

For example, let assume the filter is fed a constant input of magnitude C. What would be the best method to calculate the internal state variables as the filter approaches steady state?

When I run a simulation, it appears that S1 and S2 are negatives of one another (S1 = -S2). Is there a way to calculate the exact values?

My best guess is that I need to set-up a system of equations and then find the steady-state output (given constant input C). However, my calculations seem to fall apart when I actually try to doing them.

Edit: So after working this problem a little longer, I was able to determine the following:

$$
S_{1}[n] = X[n]\cdot (Gain_{DC} - b_{0})
$$

The DC gain is then found by taking the Z transform of difference equations and setting Z = 1.

Consider a discrete-time LTI system (that's your second order stage) described by a LCCDE in the form of a recursion equation:

$$y[n] = a y[n-1] + b y[n-2] + c x[n] + d x[n-1] $$

where $a,b,c,d$ are real or complex constants (coefficients of your IIR filter) and the values $y[n-1],y[n-2]$ are previous values of the output, where as $y[n]$ is the current value of it. And $x[n],x[n-1]$ are the current and previous values of the input respectively.

In your case, you assume a constant input $x[n] = C$ for all $n \ge 0$, hence you have $x[n] = C$ and $x[n-1] = C$ as $n \to \infty$. Now assuming that the filter has reached its steady-state output implies that $$y[n] = y[n-1] = y[n-2] = K$$ as $n \to \infty$

Hence putting those values into the recursion equation yeilds a relation between output $K$ and input $C$ as:

$$K = a*K + b*K + c*C + d*C$$ which can be arranged to yield:

$$K(1-a-b) = C(c+d)$$

At this point you tell us that output is $K=20$ mm, then you can find the corresponding input as $x = C = 20*(1-a-b)/(c+d)$, and your states will be:

$$ y[n-1] = y[n-2]= K , x[n] = x[n-1] =x[n-2] = C$$


Computing cross-correlation of two images using OpenCV I have a program that uses OpenCV to compute either the convolution or cross-correlation of an image with a specified kernel. I compute cross-correlation by setting the conjB flag to true when calling cv::mulSpectrums.

Convolution consistently produces correct results. However, the result of cross-correlation is always wrong.

Based on the convolution results, the spectra of the kernel and image are correct. The only change is that the kernel spectrum is conjugated before multiplying. Both were computed using cv::dft, and so are in the CCS-packed format expected by mulSpectrums.

I'm comparing the result with the output from scipy.signal.correlate, using Python. I tried flipping the the kernel before computing the DFT, and that produces the correct cross-correlation.

cv::dft( dft_buf, dft_buf, 1, img_in.rows );
cv::mulSpectrums( dft_buf, kernel_dft, dft_buf, 0, CORRELATE );
dft( dft_buf, dft_buf, cv::DFT_INVERSE + cv::DFT_SCALE, kernel.rows/2 + img_in.rows );

int col_offset = ( kernel.cols % 2 ) ? ( kernel.cols/2 ) : ( kernel.cols/2 - 1 ); // is odd?
int row_offset = ( kernel.rows % 2 ) ? ( kernel.rows/2 ) : ( kernel.rows/2 - 1 );

cv::Mat result_roi( dft_buf, cv::Rect( col_offset, row_offset, img_in.cols, img_in.rows ) );
result_roi.copyTo( img_out );

On closer inspection, I discovered that the erroneous correlation result resembles the correct result, but shifted up and to the left. The former was displayed in scientific format, so it was hard to see the pattern at first.

The reason is that taking the conjugate is equivalent to flipping the whole zero-padded kernel, and not just the original kernel (which is what we want).

The naive response would be to simply select another ROI from the output. However, the result is actually shifted circularly, so that wouldn't work. The correct result is present, but disjoint.

There are two possible solutions: flip the kernel before zero-padding it and computing the DFT, or change the location of the image in the zero-pad buffer.

What I did was place the image with the upper-left corner at (kernel.cols-1, kernel.rows-1). The kernel is at (0,0), but taking the conjugate flips it vertically and horizontally. Multiplying the DFTs implicitly flips it again (convolution).

Because the DFT multiplication property is equivalent to a circular convolution, the kernel buffer wraps around diagonally. So it will be the trailing portion, or "wrapped" part of the kernel buffer that "slides" across the image, and not the leading edge.

2D circular convolution can be visualized as a plain containing infinite adjacent copies of the zero-padded image. Rather than wrapping the zero-padded kernel around, it can be visualized as linear convolution with the lower-right copy, with adjacent copies providing the "history."

111100000 | 111100000
111100000 | 111100000
000000000 | 000000000
000000000 | 000000000
----------------------
111100000 | 111100000
111100000 | 111100000
000000000 | 000000000
000000000 | 000000000

For correlation with a 2x2 kernel, you could change it to this:

000000000 | 000000000
011110000 | 011110000
011110000 | 011110000
000000000 | 000000000
----------------------
000000000 | 000000000
011110000 | 011110000
011110000 | 011110000
000000000 | 000000000

You might think that the image should actually be shifted down and to the right by one more space, so that the "wrapped" kernel buffer will start with the kernel overlapping only the upper-left element of the image. The idea being that, convolution will first flip the conjugate kernel to look like this:

110000000
110000000
000000000
000000000

However, this is incorrect. The conjugate symmetry of the Fourier Transform states that flipping a signal, or an image, is equivalent to taking the complex conjugate of its Fourier Transform. However, the computer representation of the time-series is different than expected. For a 1D signal:

[ 1, 2, 3, 4, 0, 0 ]

flips to:

[ 1, 0, 0, 4, 3, 2 ]

and not:

[ 0, 0, 4, 3, 2, 1 ]

Carrying this over to 2D, our 2x2 kernel would actually look like this, after conjugating its DFT:

100000001
000000000
000000000
100000001

It's difficult to visualize how to index the correlation using this kernel. However, in 1D, this can be interpreted as flipping the array, and then circular shifting by one to the right. In 2D, it flips the kernel vertically and horizontally, and then circular shifts it down and right by one. What this means is that the kernel will be up and to the left by one more than expected at the start of convolution. I'll scale the kernel by 2 to make this easier to show. So instead of this (kernel 2s, image 1s) :

000000000
022000000
022111000
001111000

We have this (note that the image has been moved accordingly):

220000000
221110000
011110000
000000000

Remember that I'm showing the wrapped portion of the kernel buffer here, at the start of convolution.

This is hard to really show without an animation, but this is the best I could do here.

Here's the code:

// copy image to buffer
int col_offset = ( CORRELATE ) ? ( kernel.cols - 1 ) : ( 0 );
int row_offset = ( CORRELATE ) ? ( kernel.rows - 1 ) : ( 0 );
cv::Mat image_roi( dft_buf, cv::Rect( col_offset, row_offset, img_in.cols, img_in.rows ) );
img_in.copyTo( image_roi );

// Compute DFT of image, and multiply it with kernel DFT
int non_zero_rows = ( CORRELATE ) ? ( img_in.rows + kernel.rows - 1 ) : ( img_in.rows );
cv::dft( dft_buf, dft_buf, 0, non_zero_rows );
cv::mulSpectrums( dft_buf, kernel_dft, dft_buf, 0, CORRELATE );

// Compute IDFT
dft( dft_buf, dft_buf, cv::DFT_INVERSE + cv::DFT_SCALE, kernel.rows/2 + img_in.rows );

// Select ROI from result, and either copy or pass reference to output buffer
col_offset = ( kernel.cols % 2 ) ? ( kernel.cols/2 ) : ( kernel.cols/2 - 1 ); // is odd?
row_offset = ( kernel.rows % 2 ) ? ( kernel.rows/2 ) : ( kernel.rows/2 - 1 );
cv::Mat result_roi( dft_buf, cv::Rect( col_offset, row_offset, img_in.cols, img_in.rows ) );
result_roi.copyTo( img_out );


Is there an OBJECT BASED Image Segmentation method, with pure non-segmented training images? I wonder if there is a method works as follows:
1) The algorithm is trained with some images of a particular object (I train my Image with 1000 banana images)
2) I use the training data to segment a test image (semantic segmentation)

The critical point here is, the training images DO NOT contain labeling information (foreground - background). I expect the algorithm to learn the shape of the object with some common edge orientations or features etc.

I study over Graphical Models for Segmentation and the most similar method is "OBJCUT" proposed by the following paper:

http://eprints.pascal-network.org/archive/00006963/01/kumar10.pdf

But, they train the system with shape exemplars (ground truth images contain segmentations done by humans).

My objects to segment are non-articulated objects.

This is not a complete answer (I would comment directly your question but not allowed to due to rep), but you could look at "Learning and detecting shape models v1.3 (September 2009)" (http://groups.inf.ed.ac.uk/calvin/software.html). The paper's name is: "From images to shape models for object detection" by Ferrari, Jurie & Schmid

The idea is to learn a contour based model of object which is then possible to deform. Actually, you have to segment a little bit by providing simply bounding boxes (x-y coordinates of the top-left and bottom-right corner) of the object in the training images. You may also use training images that are already cropped around the object which you want to learn the model.

I tried their Matlab implementation and creating the bounding boxes for my training images was about as fast as manually cropping them in Gimp or else...

EDIT: BTW, if you wonder how to create groundtruth files with this program, I used a text editor such as "Notepad++" to create a text file. In the file, I wrote, separated by spaces, the coordinates of the bounding box, e.g.:" 10 33 572 458 " (no quote) then saved the file with an appropriate name (check how they do for their sample files) with the dummy extension ".groundtruth", their Matlab implementation will recognize the files.


Confused over Stopband attenuation and stopband ripple I have designed this bandstop filter to a desired specification:

Stopband attenuation = 40db
Passband ripple = 2db
etc ...

I understand that the stopband ripple can be given by the following:

Stopband ripple = −20log10(δs)
40 = −20log10(δs)
δs=0.01

therefore stop-band ripple δs = 0.01. where I am confused is how this value can be converted into dB so that it can be compared to the stop-band ripple seen in the MATLAB simulation (which is in dB). 40dB is the stopband attenuation, not the stopband ripple (or is it?). the simulation is showing a stopband ripple of around 32dB.... is the stopband attenuated and ripple with same ? what am I missing ?

thank you for your time

I understand that the stopband ripple can be given by the following:

  Stopband ripple = −20log10(δs) 40 = −20log10(δs) δs=0.01

Where does that come from? I have never heard the term "stopband ripple".

Stopband attenuation of 40 dB means that in the stop band any level of -40dB OR LESS is allowed. How much less than -40dB or where it hits -40dB isn't a controlled variable. -50dB, -60dB etc all are fine.


What is the sequence type if ROC is an annulus in the z-plane if:

$|z|>|\alpha|$ is causal and right-sided

and:

$|z|<|\alpha|$ is anti-causal and left-sided

then, what is the annulus/donut-ring ROC?

$|\beta| < |z| < |\alpha|$

non-casual?  finite-length?  double-sided?

if the ROC of a Z-transform is a ring on the z-plane, then the associated time domain signal is :

two sided infinite length (hence non-causal)

and

stable if the ring includes unit-circle or unstable otherwise.


Filter ECG signal with Wavelet and Python I have to filter the signal of an ECG with the wavelet method with Python. The data is in a txt file. I can create my dataframe with pandas, display that with seaborn, but can not find a way to apply the filter.

The wavelet method is imposed. I do not really know how to do it.

First I tried to understand the mathematical formula to transcribe it into a python algorithm: it's clearly beyond my abilities and my knowledge (and filling this gap would take me months, even years)

I also tried using the python pywavelets library.

The problem is that I do not know what to do with cA and cD :

cA, cD = pywt.dwt([1, 2, 3, 4], 'db1')

My research on google (in English and French) directs me towards papers of scientific research which are also outside my capacities of comprehension.

I will try to answer the question after working on this school project.

First, why filter an ECG using wavelets? I had a raw signal, full of noise. The filtering of the signal using the wavelet method makes it possible to capture spatial and temporal information very important for an unusual detection.

Here is a tutorial that inspired me, it perfectly describes the role of coefficients:: https://medium.com/@andrewtan_36013/electrocardiograms-qrs-detection-using-wavelet-analysis-a1070505efee


How to calculate time of arrival using room impulse response I need to find the timestamp of arrival of a particular signal at a receiver location within a room. The known parameters are the room impulse response h for a fixed source and receiver location, the input signal x which is a linear sine chirp.

Currently the steps that i am doing are as follows:

Convolution of input signal x and the impulse response h to generate output y
Cross-correlation of y with x, on which i am doing peak picking to determine the arrival timestamp.

Is the mentioned procedure right or am i missing something?

If you do not have a very controlled room, there might be something attenuating the direct path and a reflection might then be the highest peak. In that case you should accept also a somewhat attenuated peak if it comes earlier than the strongest peak.

Coupling to structures is not that much of an issue. For example, when a sound wave reflects at the interface of wood and air, only a fraction of about 0.0007 (see Ref. 1) of the energy of the wave will be passed to the other medium. The wave will then spend a lot of time in the first medium until much of its energy is passed to the second medium, so you can imagine that the impulse response of this coupling is not very concentrated in time and doesn't disturb the peak detection much.

References:

Acoustics of wood, course in non-destructive testing of wood, ETSI Montes, ETS Arquitectura – Universidad Politécnica de Madrid, June 2005


Filtering the EEG signals I have very limited knowledge of DSP, so I apologize if my question is trivial :)

I have an EEG signal from which I need to extract different frequency bands. For example, waves in the the frequency band of 8-12 Hz.

After noise filtering, I was planning to simply feed the data into the FFmpeg bandpass filter, setting the 'frequency' to 10 Hz, and 'width' to 2 Hz. Is this the way to do it?

I need to ask this since I'm not sure how to check if I got the correct wave in the end.

It is one way to do it, yes. You will get frequencies between 8 an 12 Hz less attenuated compared to others, but both the 8 and 12 Hz components will be attenuated 3 dB ($\frac{1}{\sqrt{2}}$ the original magnitude) and frequencies just outside the band will be attenuated only slightly more.

The way depends on how much you would like to attenuate and so on.

For a proper filter design you would like to specify a passband (between 8 and 12 in your case), a stopband (where the filter should attenuate according to your specification, there must be a transition band), maximum allowed attenuation in the passband and minimum required attenuation in the stopband. The exact numbers can be hard to determine, but once you have some you can rather straightforwardly input them into any appropriate filter design method.

However, you should obtain some filtering using the FFmpeg approach, although there may be artifacts which you did not expect/want.


Find symbol duration from symbol rate How can one find symbol duration if bit rate is 100 kbs (as BPSK is used, symbol rate = bit rate)?

The symbol duration $T_p$ is related to the symbol rate $R_s$ by the simple relationship $$T_p = \frac{1}{R_p}$$
In the case of BPSK, $R_p$ is equal to the bit rate $R_b$, and $$T_p = \frac{1}{R_b}$$
For constellations with $M$ elements, each symbol encodes $k = \log_2(M)$ bits, the symbol rate is $R_p = R_b / k$, and the symbol duration is $$T_p = \frac{1}{R_p} = \frac{k}{R_b}$$


What is the effect of padding zeros to FFT on time domain ? (opposite of what is usually done) As above. I know that zero padding in time increasaes frequency resolution but I'm having a hard time to understand the opposite. How would you interpret that?

The same effect will happen on the time-domain sequence too, due to the fact that DFT time and frequency domains are exact duals apart from a scaling factor and reversal.

Therefore, by properly zero padding (to the center of) the FFT data, the corresponding time-domain signal will be interpolated; i.e., more samples from the same signal will be obtained.


Interpret group delay of a filter I keep reading that a constant group delay over a frequency band is desirable for a filter. But a constant group delay means a linear decrease in the phase delay, which would distort the signal.

Let's assume that I have a group delay of 5 [time units] over a frequency band from 100 Hz to 110 Hz.
That means that the phase delay is $-5\omega + c$.

So:

the 100 Hz component will be delayed by 500 + c rad
the 101 Hz component will be delayed by 505 + c rad
the 102 Hz component will be delayed by 510 + c rad
etc

Considering that different frequencies components have different phase delays, the signal will be distorted.

So why is a flat group delay response desirable for a filter?

This is actually pretty simple. Take a pure-tone sine wave, say

$$s(t) = \cos\left(\omega t + \phi \right).$$

Now, fix a time $t_0$ and find the phase $\psi(t_0)$ of $s(t_0)$. We get

$$\psi(t_0) = \omega t_0 + \phi$$

Now let's vary the frequency $\omega$ and keep the time fixed. The above expression becomes a linear function in $\omega$! Even though we haven't actually "delayed" the wave.

Imagine it visually. You have a point in time (not 0, that will always stay the same) and a sine wave which is, say, shrinking (frequency increases). As it shrinks, different phases of this sine wave will cross your fixed time point for different frequencies.

This is why phase response of a simple delay system is a sloped line, not a horizontal line. Its derivative, however, will be horizontal and proportional to the amount of delay we introduced in the system, i.e

$$\frac{\partial}{\partial \omega} \left( \omega t_0 + \phi \right) = t_0,$$

Which is precisely the system delay.


Region of Convergence In attached  image why does the ROC have these values for
$$
X(z) = \frac{1}{1-\frac{1}{3}z^{-1}} - \frac{1}{1-2z^{-1}}
~~~~~,~~~~~ 1/3 < |z| < 2
$$

and for
$$
Y(z) = \frac{5}{1-\frac{1}{3}z^{-1}} - \frac{5}{1- \frac{2}{3}z^{-1}}
~~~~~,~~~~~  2/3< |z|
$$

Also why is ROC for $X(z)$ between 2 values.

According to the properties of Z-transform, the ROC of sum of two signals is the intersection of their individual ROCs.

If you look at the signal
$$ x[n] = (\frac{1}{3})^n u[n] + (2)^n u[-n-1] $$

you see that it consists of two signals; one right sided $(\frac{1}{3})^n u[n]$ with a pole at $z = 1/3$ and whose Z-transform is $X(z) = 1/(1 - (1/3)z^{-1})$ and with region of convergence ROC1: $ |z| > 1/3$ and the other is a left sided one $(2)^n u[-n-1]$ with a pole at $z=2$ whose Z-transform is $-1/(1-2z^{-1})$ and with ROC2 as $ |z| < 2$. Therefore, the intersection of their ROCs is
$$ROC_x = ROC1 \cap ROC2 = \{ 1/3 < |z| \} \cap \{|z| < 2 \} = \{ 1/3 < |z| < 2 \} $$

similarly for the signal $y[n]$.


Zero Order Hold with polyphase resampler I'm trying to design a FIR filter(direct form) the coefficient of which will be used to design a polyphase resampler that does Zero-Order Hold resampling. If the direct form impulse response is triangular, then we get a linear interpolator.

Is it possible to imagine a filter response that would produce a ZOH interpolator? I have some constraints that prevent a simple MATLAB's repelem() based ZOH implementation.

I was able to achieve this with a square pulse impulse response of length I(interpolation rate). The square impulse response in relation to an equivalent sin(x)/x and linear interpolator(triangular) response is show in the figure below. This is still not a perfect design as I see discontinuities every I samples in the resampled output.

Edit: The discontinuities are addressed by changing the transition of high-low in the square impulse response, from two samples to three sample. So instead of the transitions being (high, low) in two samples, it will be (high, mean(high,low),low) over three sample and likewise for transition from low ot high.


Where to start with DSP? I have a DSP exam coming up this summer and I got an abysmal mark on an assignment earlier this year. Obviously, the lecture slides are not enough for me to understand DSP, so I am wondering where can I go to form an understanding of DSP from the ground up.
The topics are purely math-based and do not include any programming or Matlab skill
My skills in calculus are also quite bad.

so I am wondering where can I go to form an understanding of DSP from the ground up.

There is a great list of resources in this answer What Resources Are Recommended for an Introduction to Signal Processing (DSP)?  You should spend half an hour with each one and then decide which fits your own learning style the best.

My skills in calculus are also quite bad.

That's a problem. DSP is by nature fairly math heavy and if you lacking some fundamental skills, you will we miserable the entire time. It's like trying to build a house but all you have is a hacksaw and a bent hammer. Get the right tools first. At the very least make sure you are familiar with complex number and the associated math. There are plenty of resources to do that, including Khan Academy.


How is temporal pre masking possible? I was just learning about psychoacoustic audio data compression and came across some sources suggesting that a strong tonal component (masker) that occurs within a certain short timespan after a weaker component (maskee) has the effect on the brain that, it will be unable to perceive the maskee.

What I find really hard to digest is how can it affect the causality of time? Is it because the maskee signal gets picked up by the ear but gets dropped before it gets perceived in the brain due to the masker signal happening before the communication/processing delay between the ear and the brain?

Any thoughts on this?

Auditory processing has a latency of around 40ms from stimulation at the cochlea to awareness in the brain. Those 40ms appear in a number of auditory cognitive effects and mark the separation between the perceived "now" and the past. Within these 40 ms the brain can group information and produce masking effects. Pre-masking is much shorter than post-masking, because of causality, but it's not necessarily 0.


difference between mean filter and order statistic filter I am working on my paper about comparing between mean filter and ordering statistics filter, The mean filter is contraharmonic mean filter and the ordering statistic filter is alphatrimmed mean filter. So, I compare contraharmonic mean filter and alphatrimmed mean filter

and i've read a slide presentation that explain mean filter is good for removing some kind of gaussian noise (uniform noise) and order statistics filter is good for removing some kind of exponential noise and salt & pepper noise (rayleigh noise).

Can anybody here show me the book source of that statement? is it in gonzales book? Thank you so much, your answers will help me a lot
and sorry for my poor english :(

Quick Answer

The main difference of this filters is how it perform the operations.

Mean Filter

Brief Description

Mean filtering is a spatial filter, and it's a simple, intuitive and easy to implement method of smoothing images, i.e. reducing the amount of intensity variation between one pixel and the next.

How It Works

The idea of mean filtering is simply to replace each pixel value in an image with the mean (`average') value of its neighbors, including itself.
Mean filtering is usually thought of as a convolution filter and it's based around a kernel, which represents the shape and size of the neighborhood to be sampled when calculating the mean. Often a 3×3 square kernel is used like this
 $$
   K =
  \left[ {\begin{array}{cc}
                \frac{1}{9} & \frac{1}{9} & \frac{1}{9}\\
                \frac{1}{9} & \frac{1}{9} & \frac{1}{9}\\
                \frac{1}{9} & \frac{1}{9} & \frac{1}{9}\\
              \end{array} } \right]
$$

Note: Check this link for more details in mean filtering.

Order Statistics Filter

Brief Description

This type of filter is based on estimators and  is based on "order", the sense of order is about some quatities like $\operatorname{min} $ (first order statistic), $\operatorname{max}$ (largest order statistic) and etc...

How It Works

Given $N$ observations $ X_{1}, X_{2}, X_{3}, \dots X_{N} $ of a random variable $X$, the order statistics are obtained by sorting the $\{X_{i}\}$ in ascending order. This produces $\{X_{(i)}\}$ satisfying:

$$X_{(1)} \leq X_{(2)} \leq X_{(3)} \dots \leq X_{(N)}$$

where $\{X_{i}\}$ are the order statistics of the N observations.
So, an Order Statistic Filter (OSF) is a estimator $ F(X_{1}, X_{2}, X_{3}, \dots X_{N})$.

Some common filters which fit the order statistic filter framework are:

The linear average, which has coefficients:
  $$ \alpha_{i} = \frac{1}{N} $$
The median filter, which has coefficients:
$$ \alpha_{i} = \left\{ \begin{array}{ll}
     1 & i = (N+1)/2\\
     0 & \text{otherwise} \end{array} \right.$$
The trimmed mean filter, which has coefficients:
$$ \alpha_{i} = \left\{ \begin{array}{ll}
     1/M & (N - M + 1)/2 \leq i \leq (N + M + 1)/2 \\
     0 & \text{otherwise} \end{array} \right.$$

Please check the references for more details in the Todd Veldhuizen web page.


STFT calculation with Gaussian Window $$ f(t)=\exp(jat^2) \,\,\, and \,\,\, g(t)\,\,is\,\, a\,\, Gaussian\,\, Window:$$
$$ g(t)= \left (πσ^2\right)^{\frac{-1}{4}}\exp\left (\frac{-t^2}{2σ^2}  \right ) , \,\,\,\,\,\,\left \|g(t)  \right \|=1 $$
$$ $$
I want to find the STFT (Short-Time-Fourier-Transform) of f(t) and prove that:
$$ $$
$$Psf(u,\xi)=|Sf(u,\xi)|^2=\left (\frac{4πσ^2}{1+4α^2σ^4}\right )^{\frac{1}{2}}\exp\left (\frac{-σ^2(\xi-2au)^2}{1+4a^2σ^4}  \right )$$
$$ $$
I started by calculating the Fourier Tranformation of f(t) and found that $$f(t)=\exp(jat^2)\leftrightharpoons K \cdot \exp\left ( \frac{-ω^2}{4α}\right)=F(ω) $$
where K is a constant
$$ $$
I am confused with the next steps i have to follow in order to make use of F(ω) in order to calculate STFT. Do i have to use the definition of STFT?:
$$Sf(u,\xi)= \langle\,f,g_{u\xi},\rangle= \int_{-\infty}^{\infty} f(t) \cdot g(t-u) \cdot e^\left (-j \,ξ \,t  \right )dt $$
I tried to do so but i didn't manage to calculate the integral above.
Is there an easier way to calculate STFT by using any properties?Any help is much appreciated!Thanks in advance!

so MJ, my approach to using the Gaussian window for the STFT is first to use the "ordinary frequency" definition of the continuous-time Fourier Transform:

$$ \mathscr{F}\Big\{x(t)\Big\} \triangleq X(f) = \int\limits_{-\infty}^{\infty}x(t)\, e^{-j 2 \pi f t} \,\mathrm{d}t$$

$$ \mathscr{F}^{-1}\Big\{X(f)\Big\} \triangleq x(t) = \int\limits_{-\infty}^{\infty}X(f)\, e^{j 2 \pi f t} \,\mathrm{d}f$$

and start with this really cool isomorph of the Fourier Transform:

$$ \mathscr{F}\Big\{e^{-\pi t^2}\Big\} = e^{-\pi f^2} $$

and use the well known time-scaling and translation or frequency-scaling and translation theorems of the Fourier Transform to get you

$$ \mathscr{F} \Big\{ e^{a t^2 + b t + c} \Big\} = e^{A f^2 + B f + C} $$

where the constants $A$, $B$, and $C$ can be explicitly mapped from $a$, $b$, and $c$. It appears to me that the mapping is:

$$\begin{align}
A &= \frac{\pi^2}{a} \\
\\
B &= j \frac{\pi b}{a} \\
\\
C &= c - \frac{b^2}{4a} - \tfrac{1}{2}\log\left(-\frac{a}{\pi}\right) \\
\end{align}$$

and the inverse mapping (which should be self-similar) is:

$$\begin{align}
a &= \frac{\pi^2}{A} \\
\\
b &= -j \frac{\pi B}{A} \\
\\
c &= C - \frac{B^2}{4A} - \tfrac{1}{2}\log\left(-\frac{A}{\pi}\right) \\
\end{align}$$

Looks like $\Re\{a\}<0$ and $\Re\{A\}<0$ for the integrals to converge and for the $\log(\cdot)$ to be real and finite in the mapping.


Realtime fixed point IIR filter implementation in C I would like to implement a fourth order fixed point low-pass IIR filter in C (with Q15 arithmetic), and I will prepare this filter by using two stages cascaded direct form-II bi-quad filters.

I have found countless C implementations and example coefficients on the digital signal processing books, but I saw non of them is not reliable during the verification (they create overflows, does not attenuate expected frequency components, etc.).

Which C implementation is suitable for this problem? And also, is there any sample coefficients to verify this implementation?

Although this seems like a remarkably simple questions, it requires a remarkably complicated answer.
I don't think there is a "one-size" fits all solution. The best choice of algorithm will depend on what noise you can tolerate and the type of low pass (steepness & frequency). For example at 44.1 KHz sample rate a 4th order Butterworth at 10 kHz is fairly straight forward, whereas a lowpass at a 100Hz is a royal pain. In essence it depends on how close your poles are to the unit circle.
Quantization and rounding error of IIR filters are usually transfered to the output weighed by the pole-only transfer function. A 4th order Butterworth 10 kHz low pass filter has a worst case noise amplification of only 5dB, so that is not much of a problem.
However at a 100 Hz low pass (again 4th order BW) the noise gets amplified by a whopping 75 dB. If you use Q15 math, your basic noise floor is maybe at -100dB or so. After the filter, your signal to noise ratio will only be 25 dB.
That's one of the reasons why fixed point IIR filters are fairly complicated. IF you need low cutoff frequencies and half way decent signal to noise ratio, then the basic algorithms will not work. You need to look into double precision math and/or error spectrum shaping or related methods.


Maximum frequency of any mobile phone speaker Which is the maximum frequency of a speaker of any common mobile phone? I would like to use sound to repel some animals, so I need playback sampling rate to be above 44100 Hz.

The native sampling frequency of iOS devices is 44100 Hz, while Android devices can operate at either 8, 16, 24, 44.1 or 48 kHz. Support of all of these depends on Android device manufacturer, however most Android support all of them.


Ways to identify if music has vocal timbres in it? The human voice is singular because it has a specific frequency spectrum and fast transitions between vowels and consonants.

What are good and computationally fast ways to detect if there is probably a voice contained in a noisy music, the amount of the song where there is vocals, and the type of vocals in different kinds of music, that perhaps have beats and guitars and other sounds mixed into the vocal?

Scalable audio separation with light Kernel Additive Modelling

http://www.loria.fr/~aliutkus/kaml/


Avoiding "RuntimeWarning: invalid value encountered in log10" when converting FFT data to dBFS I've attached some code below that I'm having some issues with.  Basically, I need to write a function to get the dBFS value of a signal at a specified Hz value, at a given channel in a .wav file.  For a number of reasons, I need this to work on both 16 bit integer and 32 bit integer formats.  Since SciPy's .wav io module supports 16 bit and 32 bit integer, I'm hoping that I could just convert either to 32bit float first, then perform the FFT and dBFS conversions.  For example, if you have a stereo test file that has a 440 Hz tone on channel zero, the function below would be called by db_val = check_db_at_freq_channel(test_file, 440, 0) where test_file is just a file path.

I have a test file that was exported as both 16 bit int, and a copy at 32 bit int, when I run this function on the 16 bit file, it's fine, but when I run it on the 32 bit file, I get RuntimeWarning: divide by zero encountered in log10.  Obviously, there are 0's being passed to the log10 function to convert to dBFS.  I'm wondering if there is a way I can avoid this?  Why does this only happen to 32 bit int converted to 32 bit float?  Is my method of converting to 32 bit float incorrect?

import numpy
import scipy.io.wavfil
import scipy.fftpack

def check_db_at_freq_channel(input_pcm_file, hz, channel):
    data_sample_rate, data_pcm = scipy.io.wavfile.read(input_pcm_file)

    data_len_rounded = round(data_pcm.shape[0] / float(data_sample_rate))

    if data_pcm.dtype == 'int16':
        max = float(2**15)
    elif data_pcm.dtype == 'int32':
        max = float(2**31)
    else:
        raise TypeError('data_pcm.dtype {} != {} or != {}'.format(data_pcm.dtype, 'int16', 'int32'))

    # Convert Integer data to floating point.
    data_pcm = data_pcm.astype(numpy.float32, order='C') / max

    # Check for a Mono Input
    if len(data_pcm.shape) == 1:
        data_pcm = data_pcm[:, None]

    data_fft = scipy.fftpack.fft(data_pcm[:, channel], n=int(data_len_rounded * data_sample_rate))

    # Convert to dBFS
    data_db = 20 * numpy.log10(2 * numpy.abs(data_fft) / float(len(data_pcm)))

    # Get Frequency Bins
    data_freqs = scipy.fftpack.fftfreq(data_fft.size, d=(1 / float(data_sample_rate)))

    # Get Index of requested Hz value.
    data_index = numpy.where(data_freqs == data_freqs[(numpy.abs(data_freqs - hz)).argmin()])[0][0]
    return data_db[data_index]

if __name__ == '__main__':
    test_file = '/path/to/some/file.wav'
    val = check_db_at_freq_channel(test_file, 440, 0)
    print val

I'm sure there is more information I could provide, so please let me know if I'm missing some crucial information?

Mathematically speaking, log(0) is undefined. The standard behavior for most computer languages, such as matlab, is to return -inf. Some languages, such as C++ will return -inf and raise a div by zero exception. Given that it's just a warning you could opt to ignore or suppress the warning and move on. Another option is to massage all of the 0 values to a very small positive number.

fft_mag = math.abs(fft_data)
fft_mag[fft_mag == 0] = np.finfo(float).eps


Why do we need to increase sampling frequency at the transmitter? I've thinking about this for some time now and I was wondering why do we need to increase smapling rate in the transmitter? I will explain a bit more.

From the point of view of a software-defined radio, you convert bits to symbols (let's suppose I'm doing QPSK) and after pulse shaping with 4 samples per symbol, (which I know it also increase sampling frequency by 4 since it's another interpolation after all) I've seen it's quite common to do some multi-stage interpolation (or a single-stage interpolation, it doesn't matter) to increase the sampling frequency even more. Why is this beneficial? I mean, I can only think of one reason and that'd be because you have more samples per symbol, but if after you're upconverting to passband and send it to the channel, what was that for?

Then, when you get the signal at the receiver, you decimate it to decrease sampling frequency. This one, I think I understand it better, after doing synchronization (for example in symbol sync you might need more than 1 sample per symbol), you want to decrease the sampling frequency so you don't have to deal with so many samples and is less computational costly.

So my question is, why doing interpolation in the transmitter?

Correct me if I said something wrong, please. Thank you.

Three reasons to increasing the sampling rate further are

1) To relax the requirements of the post D/A conversion filtering for image rejection.

2) Increase signal SNR by spreading quantization noise for a fixed number of DAC bits across a wider frequency range.

3) Minimize passband droop in the D/A reconstruction.

Reason 1 is the most dominant one in my opinion. The complexity and cost of a low pass filter (either digital or analog, but in this case it is the analog filter) is driven by the width of the transition band. And the reason for filtering the images is usually spectral containment as limited by regulatory requirements (and given the value of spectrum this is a HUGE reason). Reason 2 is a solid consideration if your DAC is limited by the number of bits but has wider analog BW and sufficient Spurious-Free-Dynamic Range (SFDR) to allow for overssampling benefits. Reason 3 is of least impact but I mention it as it is indeed a benefit. To address reason 3 we can often optionally simply incorporate an inverse Sinc filter digitally to compensate for the passband droop from the D/A (and this can be quite simple- see how to make CIC compensation filter).

Note interpolating DAC architectures do this for you- increase the sampling rate prior to conversion to simplify the analog anti-imaging filter requirements and give you more bits effectively since the same total quantization noise power is spread across a wider frequency range, some of which is then filtered out by the anti-imaging filter.

These graphic by Walt Kester of Analog Devices in their App Note MT-017 (https://www.analog.com/media/cn/training-seminars/tutorials/MT-017.pdf) shows this very well (Thanks Walt!):

To add, there is no benefit to the performance (sideband rejection) of the pulse shaping filter by increasing the sampling rate. The pulse shaping filters rejection is completely driven by how many symbols it spans. Two samples per symbol are all that is needed. Having more samples per symbol in pulse shaping does however simplify subsequent digital filtering if a post shaping interpolator is used- but this is exactly the same benefit described above for simplifying the analog anti-imaging filter, for the same reason.


What is the formula for the frequency spectrum? A signal $f[n]$ is given, the corresponding DTFT as $F(e^{j\omega})$ and a plot of the frequency spectrum $f(t)$. Unfortunately I can't find a formula for the frequency spectrum in my documents.

When I plot $|F(e^{jt})|$ it looks exactly like the given plot of the frequency spectrum f(t) but I don't feel comfortable to "guess" a formula by comparing plots. Searching for this in the internet confuses me more than it helps.

What is the formula for the frequency spectrum?

You are likely looking for the "Power Spectral Density" which is determined by taking the complex conjugate multiplication of the Fourier Transform. The complex conjugate process will give the expected real result and represents a power quantity (which has no phase). For the DFT this is often presented in normalized fashion such as dBFS (dB relative to full scale) or dBc (dB relative to the total power in the signal).  If the time domain signal is in volts, then scale the DFT by $1/N$ and the complex conjugate multiplication will represent the power in each bin in terms of Watts assuming a 1 ohm load (or adjust further for true power if the resistance is known). Since it is now a power quantity to convert to dB use $10Log_{10}(F(\omega)F^*(\omega))$. Where $F(\omega)$ represents the DFT or DTFT result with proper scaling, depending on which units are used.

Further as a "density" the result will be power per unit of frequency (such as dBc/Hz). The DFT and DTFT without any further windowing has an equivalent noise bandwidth that is the inverse of the time duration $T$ of the time domain signal, so therefore is equal to $1/T$ (in whatever units are used for $T$). Windowing complicates this further, but see this post for more details on that:

How to calculate resolution of DFT with Hamming/Hann window?


FFT: Removed padded zeroes? I have an application where I do this:

DFT->Filter->IDFT on a range

For computational performance I'm zero-padding my FFT to a power of two, but when I take the IDFT I want to get out an interpolated image instead of the zero-padded image. A minor detail, is that the IDFT is performed over only a certain frequency window.

What can I do in the Fourier domain to stretch my signal so that the IDFT returns the image without the zeroes I padded it with?

To make it clear, you padded your image in what domain ?

Time domain >> No problem, IFFT and then truncate will bring you the original image. It's still intact.
Freq domain >> Padding in the Freq domain means image resizing. You cant get your original image with 100% anymore. If you really want the original size one, the only thing you can do is to resize it back with some resampling methods.


Cross power spectral density (CPSD) comparison between signals I want to compare various signals with Matlab's CPSD to identify shared frequency components. (If there's a better method, I would like to know!)

I have 3 signals, A, B, and C. B is pure noise while A and C share nearby frequencies.

I noticed that CPSD's output is scaled by the power input of the signal - so for comparison between A/B/C, they should have normalized power so the CPSD amplitudes are comparable to each other (hopefully?).

The issue is that if I normalize the power, the B noise signal is boosted above the true 4 Hz signal and it looks like there's minor 5 Hz stuff for the B noise signal as well - which are both false.

How can I compare A/B/C effectively and see the relations but avoid boosted noise amplitudes in CPSD output?

I've also included a periodogram which shows the stronger 5 Hz signal but not the 4 Hz signal, for reference. mscohere gave gibberish as well. Finally, I want to keep this like "exploratory" and pretend I don't know there's a low-level 4 Hz signal that I can filter for specifically.

Any suggestions would be greatly appreciated, as I could be going in a very wrong direction or doing something incorrectly and have no idea!

t = (0:6:1080)'; %sec, some time steps
rng(1863)
A = 2*sin( 1/5*t ) + 0.25*sin( 1/4*t) + 2*(rand(length(t),1)-0.5); %some data with noise  + 2*(np.random.uniform(size=t.size)-0.5)
B = 1*(rand(length(t),1)-0.5); %some data with noise
C = 0.75*sin( 1/4.8*t ) + 0.25*sin( 1/4.05*t) + 2*(rand(length(t),1)-0.5); %some data with noise

window = hamming(110); %some window
windowPad = padarray(window,size(t,1)-size(window,1),0,'post'); %padded to t size for periodogram
Fs = 1/(1); %1/sec, time delta in freq form

pwr_A = sqrt(1/length(A)*sum(A.^2)); %estimate power of signal
pwr_B = sqrt(1/length(B)*sum(B.^2)); %estimate power of signal
pwr_C = sqrt(1/length(C)*sum(C.^2)); %estimate power of signal

[Cxy_AvB,freqs_AvB] = cpsd(A,B,window,100,512,Fs);
Axy_AvB = angle(Cxy_AvB)*180/pi;
Pxy_AvB = abs(Cxy_AvB);

[Cxy_AvC,freqs_AvC] = cpsd(A,C,window,100,512,Fs);
Axy_AvC = angle(Cxy_AvC)*180/pi;
Pxy_AvC = abs(Cxy_AvC);

[Cxy_BvC,freqs_BvC] = cpsd(B,C,window,100,512,Fs);
Axy_BvC = angle(Cxy_BvC)*180/pi;
Pxy_BvC = abs(Cxy_BvC);

figure(1);
subplot(2,1,1);
plot(1./freqs_AvB,Pxy_AvB);
hold on;
plot(1./freqs_AvC,Pxy_AvC);
plot(1./freqs_BvC,Pxy_BvC);
xlim([0 10]);
xlabel('Periods (sec)');
ylabel('Arb. Power');
title('Cross power spectral density matlab');
legend('AvB(noise)','AvC','B(noise)vC', 'Location','northwest');

A = 1/pwr_A*A; %normalize to a power of 1
B = 1/pwr_B*B; %normalize to a power of 1
C = 1/pwr_C*C; %normalize to a power of 1

[Cxy_AvB,freqs_AvB] = cpsd(A,B,window,100,512,Fs);
Axy_AvB = angle(Cxy_AvB)*180/pi;
Pxy_AvB = abs(Cxy_AvB);

[Cxy_AvC,freqs_AvC] = cpsd(A,C,window,100,512,Fs);
Axy_AvC = angle(Cxy_AvC)*180/pi;
Pxy_AvC = abs(Cxy_AvC);

[Cxy_BvC,freqs_BvC] = cpsd(B,C,window,100,512,Fs);
Axy_BvC = angle(Cxy_BvC)*180/pi;
Pxy_BvC = abs(Cxy_BvC);

subplot(2,1,2);
plot(1./freqs_AvB,Pxy_AvB);
hold on;
plot(1./freqs_AvC,Pxy_AvC);
plot(1./freqs_BvC,Pxy_BvC);
xlim([0 10]);
xlabel('Periods (sec)');
ylabel('Arb. Power');
title('Normalized Power Cross power spectral density matlab');
legend('AvB(noise)','AvC','B(noise)vC', 'Location','northwest');

[pxx,f] = periodogram(A,windowPad,512,Fs);
pxx = abs(pxx);

figure(2)
plot(1./f,pxx)
hold on;

[pxx,f] = periodogram(B,windowPad,512,Fs);
pxx = abs(pxx);

plot(1./f,pxx)

[pxx,f] = periodogram(C,windowPad,512,Fs);
pxx = abs(pxx);

plot(1./f,pxx)
xlabel('Period (sec)');
ylabel('Normalized Power');
xlim([0 10]);
title('Normalized DFT of A,B,C');
legend('A','B','C', 'Location','northwest');

I've ended up at a workable solution.

To deal with the noise, average the CPSD results of many realizations of the noise.

Pxy_AvB_mat = zeros(512/2+1,nRealizations); %preallocate
Pxy_BvC_mat = zeros(512/2+1,nRealizations); %preallocate
for i=1:nRealizations
    B = 1*(rand(length(t),1)-0.5); %new realization
    [Cxy_AvB,freqs_AvB] = cpsd(A,B,window,100,512,Fs);
    Pxy_AvB_mat(:,i) = abs(Cxy_AvB);
    [Cxy_BvC,freqs_BvC] = cpsd(B,C,window,100,512,Fs);
    Pxy_BvC_mat(:,i) = abs(Cxy_BvC);
end
Pxy_AvB = mean(Pxy_AvB_mat,2);
Pxy_BvC = mean(Pxy_BvC_mat,2);

Use a spectral analysis tool like periodogram, pwelch, or roll your own with FFTs to verify if a signal is pure noise or not. Averaging multiple realizations of noise's spectral powers will yield a flat line - showing there are no major noise constituents. Same idea as the code above. A single realization can be misleading for spectral analysis or CPSD. Many-realizations-averaged-CPSD will still show some minor peaks with respect to the noise signal, but the spectral analysis in combination will help confirm if there are any actual frequency components (noise on its own is flat).


Does it make sense to have complex numbers representing real-world audio signals? In this course (Coursera: Audio Signal Processing for Music Applications) the professor uses an example of obtaining the DFT of a complex sinusoid:

\begin{align}
x_2[n]&=e^{\jmath\left(2\pi f_0 n+\varphi\right)}\quad \textrm{for}\quad n=0, \ldots, N-1\\
X_2[k]&=\sum_{n=0}^{N-1}x_2[n]e^{-\jmath2\pi kn/N}\\
&=\sum_{n=0}^{N-1} e^{\jmath\left(2\pi f_0 n+\varphi\right)}e^{-\jmath2\pi kn/N}\\
&=e^{\jmath\varphi}\sum_{n=0}^{N-1}e^{-\jmath 2\pi\left(k/N-f_0\right)n}\\
&=e^{\jmath \varphi}\frac{1-e^{\jmath 2\pi\left(k/N-f_0\right)N}}{1-e^{\jmath 2\pi\left(k/N-f_0\right)}}\,,
\end{align}

where the frequency is expressed as $f_0$ and it has an initial phase.

So it seems like real-world (discrete) audio signals might have complex values when being represented digitally, but does this make sense?
If yes, then how can we interpret this?
Isn’t the value of the input supposed to represent a change in air pressure?
So if for an input signal $x$, say $x[0]$ was $2 + 3 \jmath$, what would this tell us about the pressure on the diaphragm at that point of time?

So it seems like real-world (discrete) audio signal might have complex values when being represented digitally,

No, you misunderstood that. The discrete audio time signal doesn't have non-real values. The Fourier transform can have such.

  but does this make sense?

It doesn't need to make sense. It's just math. I sometimes remind myself of that – it helps to occasionally get some distance from the very human desire to put sense into things.

However, assume any signal has Fourier transform $S(f)$, and that transform is purely real. Now you shift that signal in time by $\Delta t$. From knowing the Fourier transform properties, or by just applying the integral (continuous FT) or the sum (discrete FT), you instantly see that this means multiplication with a $e^{j2\pi \Delta t f}$. In the general case, this produces non-zero imaginary parts where there were none before.

In other words: even if you have a signal whose FT is purely real, simply by shifting it in time a bit, you can always get a non-zero imaginary part in the spectrum. That applies to all signals, be it audio, or I/Q baseband of a radio receiver.

  If yes, then how can we interpret this?

I think it makes sense for you to just continue to follow these lectures, I don't want to spoil too much, but:

only time-symmetrical (to be exact, time-hermitian, but since audio signals are real, that's identical to symmetry) signals have a real spectrum. If you see a non-real spectrum, you can simply tell the signal wasn't hermitian (symmetric).

Also, it tells you that if you want to deal with the spectra of signals, you always have to consider both, amplitude and phase, to represent a signal. And neither of these two aspects can be defined on the real part of the spectrum alone! The first is the root of the sum of squares of the real and the imaginary part (so the imaginary part is important), and the second is the arcustangens of the ratio of imaginary and real part. Phase "contains" the time-shift, measured in full signal periods, for harmonic signals, by the way, which is a class of signals very important for Audio.

  Isn’t the value of the input supposed to represent change in air pressure?

yes. But that the time-domain signal, not the frequency-domain signal.

  So if for an input signal x, say x[0] was 2 + 3j, what would this tell us about the pressure on the diaphragm at that point of time?

You're confusing a real time-domain signal with its Fourier transform. A single scalar value, like the momentary air pressure, is real. there's no $+3j$; there can't be, in time domain. However, in frequency domain, there might be an imaginary part.

Make sure that you understand this:

A signal is just a changing value over time, or some other axis.
The most intuitive representation hence is one along that axis, for example, an audio signal might be a change in pressure over time.

However, there's other ways of expressing the same signal. The Fourier Transform of a signal describes the same signal, but is a different function than the time signal. We often speak of "frequency domain" (after FT) and "time domain".

I don't know your background, but early on in my studies, we were taught of vector spaces, bases and base transforms.

"Time steps" is just one basis with which you can describe a signal; "Fourier Coefficients" is another one. The same signal has different coefficients (ie. representations as vectors of numbers) regarding these two bases. And the discrete FT is really, really, really nothing but a plain, old, boring base transform matrix with a few nice properties!