What Is the Difference Between MRF and Total Variation in Noise Removal? I have a general question in image processing. I have a noisy image. I would like to classify the noisy image into some regions. Two famous approaches I can use are:

MRF/Gibbs MRF: Model the spatial dependence between neighborhood
pixels.
Total Variation: Key idea may be based on smallest variation
of the image.

My question is: Could you tell me what is different between two approaches for noise removal? What is the key idea of TV? Which one is better? Thanks.

These are two different concepts that you talk about. First, MRF gives you a framework to do discrete optimization of problems, which respect the Markovian property, that is a pixel is conditioned only on the neighboring ones (roughly stated). Typical applications include binary or multi-class labeling problems. Total variation on the other hand, is generally used as a regularization by adding the integral of the absolute gradient of the signal/image to the energy functional. This helps to neglect irrelevant detail and focus on important ones.

You cannot say one is better than the other, as they are not exactly contradictory things. It depends on the application and the energy function you use in the MRFs.

An Introduction to Total Variation for Image Analysis is a good total variation tutorial to begin with. Also, Chambolle provides a link between total variation and binary MRFs, leading to new algorithms in the work entitled Total Variation Minimization and a Class of Binary MRF Models:

  Total Variation Minimization and a Class of Binary MRF Models, Antonin Chambolle, 2005


LMS Algorithm: How to Define the Desired Signal I am working on LMS algorithm,i am not getting exactly how we get desired response of the filter that is compared with the estimated output

as the others have said you are best looking at practical examples of uses of LMS. once examples is noise cancellation. think of a system where we have two microphones, one mic is the source which contains speech and background noise. the other microphone will just contain noise. in this system we try to eliminate the noise from the speech. so the inputs to the LMS are our noise reference signal and the difference between the (speech + noise) - estimate of the noise, so in this case the LMS adapts to noise that appears at the source.

Another example is echo cancellation. here you using your landline, which goes through a hybrid circuit, the signal from the far end x(n) is mixed with the near end signal u(n), so the signal going to far end is x + u + v = d, where v is noise. at the far end the person does not want to hear their own voice. so at the far end the user will have access to x(n), so inputs to LMS are x(n) and e(n) where e = d - y. in this case y is our output of LMS, used to remove the x from signal.

in a comms system, as mentioned above the LMS is used to update the equaliser coeffs. in comms systems they use the training sequence (s[n]) which is known at both ends of the transmission system, this allows us to train the equaliser. this training sequence is used before speech transmision, and during speech transmission, and sometimes training is turned off completely if channel estimate is good. in this system the LMS has inputs r[n] (the recieved signal - the s[n] after passing through the channel), and e[n], which is from the s[n] - r[n], then you would update the equaliser coeffs f after doing f = f+mu*e*r.


Detecting I/Q imbalance Are there any methods of detecting, measuring, and/or estimating I/Q imbalance (from, say, an upstream quadrature mixer and/or ADC), without using or having a controlled test signal (other than whatever the RF antenna is currently picking up).

Yes!

What is done in practice is simply disconnecting (or even better, terminating properly) the RF input of the receiver.

Then, you observe empiric noise variance separately on I and Q, and normalize so that they are equal.

This assumes that NF variations are benign, but it does work. If you can, you'd of course try to increase noise power.

USRP transceiver daughterboards can be RX IQ imbalance calibrated by letting TX transmit pseudowhite noise at a slight frequency offset (to average out potential TX imbalance) and observing the intra-board crosstalk. Same would work for any sufficiently random external signal. If you have a dominant external signal, intentionally tune a bit off, and average enough.


How to simulate a continuous signal passing through a capacitor (simple coupling capacitor) the differential equation for the current flowing through a serial capacitor (see for example https://www.allaboutcircuits.com/textbook/direct-current/chpt-13/capacitors-and-calculus/) indicates that one must take a derivative of the signal but in the literature, a serial capacitor (for instance a coupling capacitor) is often considered as a unit that just eliminates a continuous component of the signal. Well. Is it necessary to take a derivative or just get the first moment of the stochastic signal and eliminate it in order to simulate the passing of that signal through the coupling capacitor?

Is it necessary to take a derivative

Yes

or just get the first moment of the stochastic signal and eliminate

That's just a very crude approximation of taking the derivative. In some cases that's sufficient, but in most it's not.
For DC blocking the capacitor is typically used together with a resistor to form a high pass filter. The approximation is often "good enough" if the corner frequency of the high pass is well below any signal frequency of interest.


Stability varying with order I have been doing some research into modelling acoustic spaces, and (with help from this site) I've managed to take the frequency response of fibrous layer model, convert it to a reflection response as a target for an IIR, and then find the polynomial coefficients for it using FDLS.  All good so far.

I performed the usual BIBO stability test of finding the magnitude of the IIR denominator roots, and with a 4th order filter everything was fine (<=1.0).  But increasing the order, even by just one, then it became unstable.

My question is: What does the order actually represent response-wise?  I always assumed it determined the sharpest change in direction of the frequency response i.e. the 'jagged-ness' or complexity of the response curve.  If that's true, then my filter is unstable, but the few poles I've calculated just happen to be inside the unit circle.

But is that true?

FDLS is similar to rational polynomial interpolation in that the interpolation quickly becomes unstable if the target polynomial degree is set too high.  (e.g. just a tiny amount of variation or noise in the spec at a sub-optimal grid spacing can cause extremely wild excursions in the graph of a high order polynomial).  Thus it is usually better if the FDLS order you choose is overdetermined by your specification points to some degree (usually selected by trial-and-error).

The jagged-ness of a frequency response is more related to how close a group of poles or zeros are to the unit circle and each other, rather than to their absolute number.


Fitting an IIR filter to a complex transfer function anybody knows of better ways to fit an IIR filter to a complex transfer function than Matlab's invfreqz.m? I'm implementing a little transfer function fitting toolbox for myself. I started with re-implementing the methods used in Matlab's invfreqz.m. Invfreqz.m is using Levy's complex curve fitting algorithm to compute starting coefficients for a iterative Gauss-Newton algorithm. Since Levy's paper is not documenting the algorithm very well and it's quite old, I'm looking for easier implementations and maybe even better algorithms to do the same job. Literature recherche turns out to be a real pain - so any hints from experienced people regarding literature, starting points or examples for implementations would be very appreciated.

As I mentioned in a comment, I think that the equation error method is a very good starting point for designing an IIR filter in the frequency domain with prescribed magnitude and phase responses. For an explanation of this method have a look at this answer and the references there.

Two problems associated with the equation error method are

stability of the designed filter is not guaranteed
the frequency domain error is weighted with the squared magnitude of the denominator polynomial of the designed filter, which means that the approximation is better in certain frequency regions than in others, but this weighting is a result of the design process and is usually undesirable.

The first problem can be tackled by adapting the total delay of the desired frequency response such that it matches with what a causal and stable filter with the chosen orders of numerator and denominator polynomials can achieve. So the desired complex frequency response $D(\omega)$ is changed to

$$\tilde{D}(\omega)=D(\omega)e^{-j\omega\tau}\tag{1}$$

with some delay parameter $\tau$ (which can also be negative). Optimizing the value of $\tau$ will help designing a stable filter with the smallest possible error. The only difference with the original desired response is an added (or subtracted) delay (i.e., a linear phase component).

One option to address the second problem is to solve a sequence of equation error problems with an adapted weight function, which is the basic idea of the Steiglitz-McBride method.

I've implemented the equation error method for the frequency domain design of IIR filters, and I've also added the possibility to iteratively approximate the optimal least squares solution, as described above. Note that convergence cannot be guaranteed, but in practice convergence is usually achieved after a few iterations. Also note that the resulting filter may be unstable, but this problem can be solved by adjusting the total delay of the desired response, as explained above. The Matlab/Octave code can be found at GitHub.

Here is a simple filter design example. The desired response is a low-pass filter with a linear phase response in the pass band. The delay is desired to be an integer number of samples, but the exact value is to be optimized to get the best approximation by a stable filter. The numerator and denominator orders are chosen to be $M=N=10$. Experiments showed that the best value for the desired pass band group delay is $28$ samples:

w = pi*[linspace(0,.2,20),linspace(.25,1,75)];   % frequency grid
D = [exp(-j*w(1:20)*28),zeros(1,75)];           % desired response
W = [ones(1,20),10*ones(1,75)];                  % error weighting
[a,b] = eqnerror(10,10,w,D,W,6);                 % 6 iterations of equation error method

The figure below shows the design result. Note that the pass band group delay approximates the desired $28$ samples. The designed filter is stable and has a maximum pole radius of $0.98$.


Additive disturbance in block diagram representation Drawing the open loop block diagram between valve and pump is straightforward (please see illustration). But why should $v$ be an additional disturbance in this system (please see solution)?
$v$ is going out of the tank?

The pump keeps the tubing pressurized, the desired input flow is tuned via the valve. However the net flow rate that feeds the tank is not this inflow, it is the difference between the inflow and outflow.
The outflow caused by the pump is modeled as a disturbance because, in the current example, the goal is to model the liquid's level inside the tank. The outflow maybe considered a byproduct of the pump's operation.
Although the above diagram describes the plant (the physical process of the hydraulic system), if the control of the plant was studied later a possible design goal could be to maintain a constant and controllable level of the liquid in the tank (this requires feedback control aka a closed loop) and/or reaching the steady state value(desirable level in the tank) fast.
In addition a pump's flow rate tends to decrease as the pressure inside the tubing increases.
Think about a garden's hose, which sinks water into a flower pot at some flow rate and the water drains out of the flower pot through the bottom. Then if the inflow is much larger than the outflow the soil's surface will be flooded not being able to properly absorb the incoming flow of water. On the other hand, if the outflow is close enough to the inflow then the soil will be able to absorb the incoming water flow, not resulting in the flower pot's water overflow.
The block $G_t(s)$ describes the tank's response therefore the inflow-outflow difference is its input. The system $G_v(s)$ is the valve's response to the user input.


calculate fft convolution for an LTI system given its input (in time domain) and its poles and zeros I have considered taking the fft (MatLab) of the input (x) and create frequency bins for the first half (problems with an odd number of samples, tried solution zero paddings) and evaluate the transfer function in these freq bins. Then apply conjugated symmetry for completing the fft of the Transfer function. This pair of sequences are multiplied term by term and then apply ifft (MatLab). Finally scaled by the sum of the fft signal of the transfer function.
There is something strange/invalid with this procedure? How I can correct it? Bibliographical references are appreciated.
Thanks, again.

In general the filter you want to implement is an infinite impulse response (IIR) filter, unless all poles are at the origin of the complex plane (assuming causality), in which case it is a finite impulse response (FIR) filter.
In the general (IIR) case, your suggested method will not result in an exact implementation of the filter. There are two reasons for this:

by using an FFT of the filter's frequency response you actually approximate the infinitely long impulse response by an impulse response of finite length. That impulse response is simply given by the inverse FFT of the sampled frequency response.

multiplication of the FFTs implements circular (cyclic) convolution, which is different from linear convolution.

Both errors can be made small. For the first error (FIR approximation of IIR filter), you just need to choose an FFT length that captures a large percentage of the energy of the impulse response. This basically means that you choose an FIR filter of sufficient length to approximate the given IIR filter. For minimizing the second error (circular convolution instead of linear convolution), you need to zero-pad the input sequence and the FIR approximation of the IIR filter.
The question remains why one would want to use an implementation like this. One disadvantage of the proposed solution is that you have to wait for the complete input signal before you can start processing, i.e., you introduce a substantial delay. This problem could be tackled by block processing, like overlap-save or overlap-add in the case of standard fast convolution FIR filtering. The other disadvantage is the increase in memory requirements and computational load. Most practical IIR filters have relatively low orders (lower than $20$), but an FIR filter that provides a reasonable approximation will usually have many hundreds of coefficients, or even more.
There exists an exact method for block processing of IIR filters, in which FFTs can be used to solve certain matrix-vector multiplications. This method is explained in detail in this very accessible book chapter by Selesnick and Burrus: Fast Convolution and Filtering (Section 8.3.1).


Fast phase calculation I need to detect the PSK modulation posted on my previous post(phase difference detection) with 8 to 16 phase constellations (depending on selected baud rate) over a max 60 degree PSK phase range.For communication start and end the Tx generates carrier with delta(ph)=0 or -180 deg to enable Rx to synchronize or detect end of communication.as speed of converging to the result is important to my application (theoretically real-time application) I was wondering if there is a fast algorithm with non- or low-iterative mechanism that can converge fast to the instantaneous and differential phase value (implementation with verilog HDL).

I think the sin(delta(ph)) value can be filtered out of I/Q component multiplication relatively fast, but it needs an LUT to save the sin(delta(ph)) values for all phase constellations. new ideas are always welcomed!

Below summarizes efficient phase estimators for this application updated to include both a phase range of +/- 30 degrees and +/- 60 degrees. This is given in two parts, estimators for a real IF (intermediate frequency) signal, and estimators for a baseband complex signal. At the end are additional considerations related to acquisition.

Efficient Phase Estimators for Real IF Signals

Product detector: For real signals, a common phase estimator (detector) is a multiplier followed by a low pass. For this application where sensitivity is desired over $\pm 30$ degrees, the signals are nominally 90 degrees in phase resulting is an estimate that is proportional to the sine of the phase between two signals:

$$y(t,\phi) =A_1\cos(\omega_ct)A_2\sin(\omega_ct+\phi) = \frac{A_1 A_2\sin(\phi) + A_1A_2\sin(\omega_ct+\phi)}{2} $$

Where when followed with a low pass filter removes the time varying component resulting in:

$$y(\phi) = \text{LPF}[y(t,\phi)] =\frac{A_1A_2}{2}\sin(\phi) $$

Where $\text{LPF}[\cdot]$ is the time average provided by a low pass filter.

As a demodulator, this can be implemented either in a coherent receiver where $A_1\cos(\omega_ct)$ is estimated during acquisition and provided as an NCO in a digital implementation (or VCO in an analog), or in a non-coherent receiver (where the interest is in the phase difference between two successive symbols) where the demodulation is done by multiplying the received signal with a time delayed copy of itself, delayed by one symbol duration plus the duration of a quarter cycle of the carrier to convert $\cos$ to $\sin$ (when the IF carrier is sufficiently larger than the symbol rate):

$$y(\phi_2-\phi_1) = \text{LPF}\bigg[A\cos(\omega_ct+\phi_1)A\cos(\omega_c(t-T_s-T_c)+\phi_2)\bigg]$$

Where $T_s$ is the symbol duration in seconds and $T_c = 1/(4f_c)$ is a quarter cycle of the IF carrier in seconds with the IF carrier frequency as $f_c$ in Hz. Resulting in:

$$y(\phi_2-\phi_1) = \text{LPF}\bigg[A\cos(\omega_ct+\phi_1)A\sin(\omega_c(t-T_s)+\phi_2)\bigg]$$

$$ y(\phi_2-\phi_1) = \frac{A^2}{2}\sin(\phi_2-\phi_1)$$
$$ y(\Delta\phi) = \frac{A^2}{2}\sin(\Delta\phi)$$

A very efficient way to implementing either of the above approaches digitally is by hard limiting the input signal which reduces the above to a simple XOR of the most significant bit of the waveform. For accuracy, this requires ensuring the inputs to the XOR operation are a 50% duty cycle, but the result is linearly proportional to phase! It is usable over a $\pm 90°$ range with a linear phase result. Further, hard limiting a phase modulated waveform provides a 3 dB SNR improvement in positive SNR conditions (since all AM noise is removed), but can be more susceptible to jamming an interference (3 dB loss in negative SNR conditions). This is an approach to be considered due to the simplicity and high phase linearity.

As above, the X-or phase detector could be used in a coherent receiver where the NCO is also simplified to a 1 bit output (basically the MSB of a counter and you increment the count rate to adjust the frequency as part of a carrier tracking loop), or non-coherently where the MSB of the received signal is XOR'd with a delayed copy. As with the multiplier the two input signals would be in quadrature to center the detector over its unambiguous range.

Efficient Phase Estimators for Complex Baseband Signals

Given a generic signal as $Ae^{j\phi}= I + jQ$, the actual phase is given by $\phi =\tan^{-1}(Q/I)$ or $\phi =\sin^{-1}(Q/A)$.
The following is a summary of efficient phase demodulation approximations for a variation over $\pm30$ degrees and $\pm60$ degrees, assuming carrier recovery and timing is established during the $0 / 180$ acquisition period. Initial thoughts on approaches to efficient acquisition are also included at the bottom of this post.

Summary of Results

Below is a table summarizing the peak and rms phase error for various estimators. Estimators that were included in earlier versions of this post that offer no advantage to those listed below have been removed. As Ben suggests in the comments, the Q/A estimators are attractive for FPGA implementation since A is assumed to be constant over the duration of the packet.

Plots showing the relative performance are included below:

Detailed Descriptions

The estimators that are scaled by the envelope magnitude $A$ (Q/A, Q/A Juha and Q/Est(A)) are preferred since $A$ can be readily determined during acquisition of the 0/180 signal, and only needs to be determined once for relatively short packets, or is a parameter from the AGC otherwise. In a constant envelope phase modulated signal such as this, the received signal can be simply hard limited if there isn't a concern with potential 3dB loss with stronger out of band interference (or the complete loss from hard-limiting in the presence of a coherent jammer). Further, there is no need to actually divide by $A$, assuming $A$ is maintained to be constant over the packet duration the result will be linearly proportional to the phase and the decision thresholds can be set accordingly.

Q/A

$$\phi =\sin^{-1}\bigg(\frac{Q}{A}\bigg)$$
$$\frac{Q}{A} = sin(\phi)$$

for small $\phi$,  $sin(\phi) \approx \phi$ for $\phi$ in radians:
$$\phi \approx \frac{Q}{A}$$

Q/A Juha

Similar to  @JuhaP's suggestion in the comments of removing the linear slope error for the $Q/I$ estimator, here applied to the Q/A Estimator. The coefficient is found from the linear portion of the remaining terms in the Taylor Series expansion that weren't used, minimizing the error:

For ±30° Operation:
$$\phi \approx 1.0475\frac{Q}{A}$$

For ±60° Operation:
$$\phi \approx 1.150\frac{Q}{A}$$

Q/Est(A)

A fast and very efficient approach for estimating magnitude is the $\alpha$ max plus $\beta$ min algorithm where the maximum between $|I|$ and $|Q|$ scaled by coefficient $\alpha$ is added to the minimum scaled by coefficient $\beta$. At 30° range, $Q$ would always be the minimum and  $I$ always positive so this would simplify to $\alpha I + \beta|Q|$. A common choice for FPGA implementation is $\alpha = 1$ and $\beta =1/2$ since this minimized the error over all phases with bit shift divisions, but in this case $\alpha = 1$ and $\beta =1/4$ is a better choice given the narrowed phase range of $±30°$. If  multipliers were acceptable, the optimized coefficients are $\alpha = 0.961$ and $\beta =0.239$.  The plot below summarizes the two choices:

$$\phi \approx \frac{Q}{\alpha I + \beta |Q|}$$

option 1:  $\alpha =1$, $\beta = 0.25$
option 2:  $\alpha =0.961$, $\beta = 0.239$

Also not plotted below but shown above is the option optimized for use over ±60°:
 $\alpha =0.85$, $\beta = 0.45$

Note these are not optimized for the estimate of $A$, but to minimize the phase estimation error.

Q/I Phase Approximation

$$\phi =\tan^{-1}\bigg(\frac{Q}{I}\bigg)$$
$$\frac{Q}{I} = tan(\phi)$$

for small $\phi$,  $tan(\phi) \approx \phi$ for $\phi$ in radians:
$$\phi \approx \frac{Q}{I}$$

As @JuhaP mentioned in the comment, the linear slope component of the error could be removed by multiplying by 0.9289 resulting in (This one is labeld Q/A JuhaP in the plot). The coefficient below is slightly different than his suggestion but minimizes the error as it was found from the linear portion of the remaining terms in the Taylor Series expansion that weren't used rather than his approach of a first order terms of a polynomial fit to arctan:

$$\phi \approx 0.9289\frac{Q}{I}$$

Taylor Series Phase Approximations

The first term is the Q/A and Q/I approximations covered above for $\sin^{-1}$ and $\tan^{-1}$ respectively. Going beyond that is NOT recommended if efficiency is paramount but included for accuracy comparison.

arcsin

$$sin^{-1}(n) = \sum_{n=0}^\infty \frac{2n!}{2^{2n}(n!)^2}\frac{x^{2n+1}}{2n+1} \text{ for } |n|\le1$$
$$sin^{-1}\bigg(\frac{Q}{A}\bigg)= \frac{Q}{A} +\frac{1}{6} \bigg(\frac{Q}{A}\bigg)^3 +\frac{3}{40}\bigg(\frac{Q}{A}\bigg)^5 ... \text{ for } |Q/A|\le1$$

arctan

$$tan^{-1}(n) = \sum_{n=0}^\infty (-1)^n\frac{x^{2n+1}}{2n+1} \text{ for } |n|\le1$$

$$tan^{-1}\bigg(\frac{Q}{I}\bigg) = \frac{Q}{I} -\frac{1}{3} \bigg(\frac{Q}{I}\bigg)^3 +\frac{1}{5}\bigg(\frac{Q}{I}\bigg)^5 ... \text{ for } |Q/I|\le1$$

Using the first two terms for each results in:
$$\phi \approx \frac{Q}{A} +\frac{1}{6} \bigg(\frac{Q}{A}\bigg)^3$$.

$$\phi \approx \frac{Q}{I} -\frac{1}{3} \bigg(\frac{Q}{I}\bigg)^3$$.

A linear slope can also be removed from either of these with a gain constant multiplication as was done with the Q/A and Q/I estimators to further minimize the error.

Other Estimators

Juha Squared

@JuhaP offered this interesting estimator in the comments. Not very efficient but highly accurate with square terms:

$$\phi \approx \frac{3QI}{Q^2 + 3I^2}$$.

Acquisition

Efficient Acquisition for 0/180 preamble:

One idea that comes to mind for acquisition during the $0/180$ transitions is to use $\text{sign}(I_2)Q_1-\text{sign}(I_1)Q_2$ to get the change in phase between two symbols, which can the be corrected in a fast coverging and simple loop by derotating the incoming signal.  This appraoch would work well if the frequency offset is such that the phase does not rotate more than $\pm \pi/2$ between successive signals, otherwise a course FLL can be used first to get the offset within this acquisition range.
For a coherent receiver approach a PLL would be used to lock/track an NCO or PLL to the carrier and my squaring the received signal a reference tone at twice the carrier can be tracked for all modulations presented here (both the bi-phase acquisition interval and the 30 degree phase modulation when doubled will produce a distinct tone at 2x the carrier). Similarly a Costas Loop would track both signals while providing the reference signal that is nominally 90 degrees in phase with the carrier, thus providing both carrier recovery and phase demodulation.

Sources:
Taylor Series Expansions for arcsin and arctan:
https://proofwiki.org/wiki/Book:Murray_R._Spiegel/Mathematical_Handbook_of_Formulas_and_Tables


What is the bandwidth of a (real) sinusoidal tone, and pulse? I would like to know how to go about calculating the bandwidth of:

A constant (real) sinusoidal tone
A (real) sinusoidal pulse.

The question is as simple as that, but I am having a hard time with the concept of what exactly the bandwidth of a constant tone should be to begin with, and from there, what the bandwidth of a pulse should be.

In the frequency domain, a constant real tone of frequency $f$ exists as two delta functions, located at $f$ and $-f$, but how does one go about calculating its bandwidth?
Furthermore, in regards to the pulse, this is rectangular function in time, and thus a sinc in frequency domain, so would not its bandwidth simply be $\frac{1}{T}$, where $T$ is the duration of the pulse?

The spectrum of a continuous tone is, as you said, of the form $\delta(f-f_0) + \delta(f+f_0)$: 2 impulses at frequencies $f_0$ and $-f_0$.

As a lowpass signal, this is said to have bandwidth $f_0$ (the one-sided spectrum has components up to $f_0$).

As a bandpass signal, it has zero bandwidth (there's nothing around the carrier frequency $f_0$).

If you multiply the sine wave by a pulse, this makes it time-limited, and therefore frequency-unlimited. Infinite bandwidth in theory.

In practice, you must define some criteria for estimating your bandwidth. Examples are:

3 dB drop (of the sinc function around $f_0$)
10 dB drop
drop below the noise level


Can I combine an IIR (an EWMA) with a maximum? I'd like to have a exponentially weighted moving average (EWMA) which is raised to the value of the input signal whenever this is higher than the filter output.  This taking of the maximum shall take place during the iteration, not afterwards (which would be simple).

So, the algorithm should look like this:

o = 0  # or some arbitrary initial value
for i in input:
    o = (o * 99 + i) / 100
    if i > o:
      o = i
    print o

So effectively, my output shall fall slowly after a quickly falling edge in the input without rising quickly before a rising edge.

In fact, any other algorithm which just copies the input except for the places where it has a rapidly falling edge (there it should play parachute and let the output fall slower) would do fine.

The only restriction I have is that this needs to be implemented in the usual Python libraries (numpy, scipy, pandas); iterating myself will not be fast enough.

I tried applying a rolling maximum, followed by a rolling mean with a Gaussian window of the same width (using pandas.Series(arr).rolling()); this gives me nice smoothing of the fall after a falling edge; unfortunately it also gives me an equal smoothing before the rising edges (which I don't want).

If there was a way to apply an asymmetric Gaussian window (e. g. rising slowly and falling quickly) that would be a solution I guess.  But I haven't found a way to achieve this yet.

What you are describing is an envelope detector which can be formed by a Peak Detector followed by an ADSR.

The peak detection is performed by the non linear 'max' filter as per your description. This, on its own, would track the maximum of the signal and stay at that max value.

To now modify that signal's dynamics so that it 'launches' and 'parachutes' at specific rates, you add the Attack-Decay-Sustain-Release component. You don't need all four stages. This could simply be an Attack-Decay (as per your needs).

In fact, in your application, the Attack is zero. As soon as a new max is detected the output of the system is reset immediately. But, while the new values are lower than the maximum, this 'output' now decreases at a desired rate (The Decay stage). In your case this could be as simple as o-=aRate if o>0 else 0 in the else branch of the if.

Hope this helps.


Find highest frequency of a very manky signal I'm trying to compute the highest frequency (as can be sampled) in some pretty manky looking discrete time-dependent signals. My current method - a discrete fourier analysis - fails for some pretty awful looking but clearly oscillating signals (with discernable highest frequencies).

My current method is to compute the discrete fourier transform, locate the local maxima, and perform weighted averages around these peaks (to undo the 'smearing' of the dft across the discrete sample frequency bins).
Here it is in MATLAB.

% compute spectrum
dft = abs(fft(weight));         % real response
dft = dft(1:floor(end/2) + 1);  % keep only pos freqs
dft = dft/sum(dft);             % normalise

% find present modes modes
mode = zeros(0, 2);

[pks, locs, widths, ~] = findpeaks(dft);
for k=1:length(locs)

    % only consider significant peaks
    if pks(k) > 0.1

        av_width = min([floor(peak_width_factor*widths(k)), locs(k)-1, length(dft)-locs(k)]);
        inds     = (locs(k)-av_width):(locs(k)+av_width);
        av_ind   = (inds * dft(inds)) / sum(dft(inds));
        av_freq  = (av_ind - 1)/T;

        mode = [mode;
                av_freq, pks(k)];
    end
end

This produces a matrix of present modes (of a spectral significance above 0.1/1) where each row is the frequency and significance of the mode.

This works great for signals like this (time signal left, dft right with the frequencies of detected modes labeled):

but fails for signals like these (where we see and expect frequencies close to 2; not an order of magnitude smaller!)

Often I expect a certain frequency (e.g f=2.6 in the above signals) and can judge myself the mode is present by the average period between the local maxima. I tried codifying this - computing the average time between local maxima - in MATLAB, but it was pretty unreliable:

[~, peak_inds] = findpeaks(weight);
cycle_periods = diff(peak_inds) * dt;
av_cycle_period = sum(cycle_periods)/length(cycle_periods);
av_freq = 1/av_cycle_period;

My signals are well enough sampled (around 15 to 20 values per observable manky period) to resolve these manky modes visually.
I've studied time series and random processes at an undergraduate level, but we never really went too deep into spectral analysis. So:

How can I reliably compute these highest-frequency modes in my very 'unfourier' signals?
Why does my current DFT analysis incorrectly deduce very low frequency modes in these manky signals?
Why would, if multiple present modes have dissimilar frequencies, my naive distance-between-local-maxima-average method fail at extracting the highest frequency?

I think most of your questions can be solved by subtracting then mean from the signal. Namely the all sinusoidal waves, with a nonzero frequency, have a mean of zero. So when the mean of a signal is not close to zero, then it will show up at 0 Hz (you can look at this as $\cos(0\,t)=1$).

After a closer look and some testing by myself it does seem that only removing the mean would not fix your problem. Your signal seems to be contaminated by white noise, which has been integrated probably twice (so low frequencies are much more present than high frequencies). You could try to subtract a higher order polynomial fit, so linear or even quadratic. These fits can easily be obtained using some linear algebra, for example for a quadratic fit:

X = [ones(N, 1), t', t'.^2];
b = (X' * X) \ X' * y';
z = y - (X * b)';

Or you could also look at the biggest peaks above $f^{-2}$, with $f$ the frequency vector. However if you have a signal with normal white noise you should probably just stick with removing the mean, since my proposed methods try to make use of the fact that the present noise is white noise which has been integrated twice.


Sound synthesis effects for natural sounds I´m currently working on an accordion synthesizer and I would like to make it sound as natural as possible.

I really like how the following wave sounds:

https://dl.dropbox.com/u/20437903/onda%20acordeon.wav

The wave looks like this:

Looking at the wave it looks like it had harmonics added and substracted. Am I wrong?

I´m a bit new to sound synthesis and would be very grateful if anyone could please help me identifying the effects applied to the wave and some guidance to apply them. Thanks.

First, grab Sonic Visualizer, it is much better than Audacity for looking at sounds.

What you see here is probably the result of the sum of two simple and stationary sounds at fundamental frequencies close to each other. This causes beating of their fundamental, causing the amplitude modulation (tremolo) you observe.

Two important factors make a synthetic sound "feel" natural:

Whether it accurately reproduce the sound of the original instrument over time. Some criteria to look at: timbre (distribution of harmonics), timbre modulations over time, amplitude modulation over time (presence of tremolo, overall amplitude envelope), pitch modulations (presence of vibrato, does the instrument instantly reaches the target note or is there a short transition like on brass instruments...).
How it responds to control by the performer. You have to list which parameters, along with note data, will be available to your system (will it synthesize sounds from a bare MIDI file? From a MIDI keyboard sensitive to velocity and aftertouch? From a dedicated hardware device emulating an accordion, with accelerometers and pressure sensors?), and understand how to map these inputs to synthesis parameters. You also have to list which controls the performer will have over the system (for example, is there an option for the user to switch on/off the different registers?).

I will mention 3 synthesis approaches you can follow for this problem:

Sampling. Record, one by one, all the notes of an accordion, under different performance parameters (how hard the below is moved for example, all different combinations of registers). Play back from this exhaustive database of sounds, with loops to make them last as long as you want. Advantages: a note taken in isolation sounds indistinguishable from the real thing. Drawbacks: might sound "fake", and will not respond to performer input (not "expressive") unless you start doing some signal manipulation on the recorded sample - at which point they won't sound as natural as they were recorded. A lot of disk/memory used by the samples. Warning: it's quite an involved job to exhaustively sample an instrument and organize all the data, and there are already companies very good at doing this kind of stuff - it'll be hard to beat them.
Classic subtractive synthesis. You could start with two square or pulse waves, slightly detuned to get beating, a tad of high-pass filtering to give it a "nasal" quality, and an EQ to color this... Then stack several of those to get different registers. Check in the patch libraries of classic "old" synths (JX8p, D50, M1)... they all have decent accordion sounds which are totally sample-free. Advantages: can be made very expressive - since you have control over the whole production of the sound, you can easily map input parameters on sound parameters. You won't have to reinvent the wheel and can base yourself on existing sound synthesis code. Inconvenient: will sound "synthetic" (but "lively & synthetic" is sometimes better than "realistic & dead"!).
Physical modeling. You could get something more realistic than vanilla subtractive synthesis by using physical modeling - try cobbling together (for example in csound, Max/MSP or supercollider) a few reed models (the 3 mentioned programs have basic primitive for physical modeling - but I'm not sure they have free reeds), add EQ or a simple resonator. Ideally, you could get something both realistic sounding and controllable, though it's not the easiest path!


How does sampling a frequency result in a duplicated spectrum? This is probably a simple question, but I have a hard time understanding how c. converts to d. in the following picture.

Further reading as an answer is acceptable. After all, I am only asking this because the PDF I’m reading doesn't really explain why c. results in d. The more I look at it, the more I am confused.

if you consider the spectrum of analog signal as the figure in your picture.Sampling in time domain is same as multiplying a pulse train to  the analog signal.This means that the frequency response of pulse train is convoluted in frequency domain. The frequency response of pulse train is again a pulse train,hence the spectrum gets duplicated in frequency domain due to its convolution with the pulse train.


Computing the gains of a piecewise affine system I have a piecewise affine (and memoryless) system, e.g., a response like the following one.

I could manually find the gain in each region by applying a ramp signal. However, the system's input must be a sine wave.

Is there any technique that could be used to find the individual gain and slopes in each region?

Add a DC signal to your sinusoidal input, so that the input falls in the middle of each line segment of the input-output curve. Then choose the amplitude of the input sinusoid in order to avoid "spilling over" into the other segments. Measure the amplitude of the output sinusoid (neglecting its DC component, of course). Lastly, compute the gain.


How to detect if the given signal is Data or Clock signal type? I have the edge information of the signal. I need to find if it's a Clock or Data type signal without iterating through all the points in edge vector...
Please provide your valuable inputs

Clocks are very regular, i.e. all the up and down transitions are equally spaced. Data tends  to more irregular, i.e. spacing between transitions varies.


Kalman Filtering with Unknown State Transition Matrix I'm currently studying the use of Kalman filters for estimating linear systems. My current State Transition Matrix (STM) is the identity since so far I've been dealing with non time-varying systems. However my next step is to extend to time-varying systems and in this case there is no information on how the system varies and so the STM is unknown. And as of my testing so far a wrong STM causes a ton of errors on the Kalman Filter estimations. So my question is if there is a method for circumventing this, be it by updating the STM on the fly or some kind of way of updating the Process Noise Covariance Matrix in order to account for that?

Unfortunately, in Classical Kalman filter applications you need to have a guess of your state transition model governing the states. Classical choices for such models for time-varying systems are the random walk model, or its generalization, a first-order Markov model. If do not have a clue about which state model you can assume, you can try to extract a trend component from your state time series (via EMD or a moving average filter, for example), and plot it to have an idea of how the first-order moment of your states vary with time, for example. If it decays as a power law with $n$ (i.e. $E[x(n)] \approx \alpha^n$) until a final value, it may be an indication that a first-order Markov model with pole equal to $\alpha$ can serve as a state transition model.

Otherwise, you can treat the transition model for your states as a black box and consider techniques such as the Unscented Kalman Filter, for example.


Fast cross correlation with limited range I am measuring the phase shift between two pseudorandom signals using correlation ifft(A.conjugate()*B) and picking out the maximum).
That works very well, but takes a lot of time. I roughly know where the peaks are, and don't need the correlation for other shifts.

Is there an adaptive cross-correllation algorithm that would speed up the process if I am only interested in the positions of maxima, and already roughly know their positions?
Analog to the sliding DFT?

Maybe running the cross-correlation on a decimated version of data, finding the peaks, interpolating, and then doing a manual correlation around the max?

I already know within 200 samples where the maxima are.
The record length is about 200k samples.

If the limited range is significantly more than log(N) shorter than the total signal of length N, then ordinary cross correlation over that limited range may be faster than FFT-based fast correlation.  You could just iterate away from the estimated/predicted lag until you find you have gone past a correlation peak (by a sufficient amount to account for expected noise).  But if you go more than log(N)/2 points away from your initial estimate, then perhaps switch back to fast convolution cross correlation.


What Does It Mean by Gradient / Divergence / Curl of a 2D Image $ g \left( x, y \right) $? Suppose, a 2D image is represented by g(x,y).

Then, what does it mean by "Gradient of g(x,y)", "Divergence of g(x,y)" and "Curl of g(x,y)"?

Please, avoid formal math for explanations.

Please, use layman's terms and informal discussions.

Original Image

Image Gradient

The gradient shows how quickly the color (or greyscale) changes from pixel to pixel. In very "busy" images, the gradient will change much. In areas of less color change, the gradient will be close to zero.

In the image above, the petals don't change color very much so the gradient is close to zero (or small). At the edges of the petals where the color changes from purple to green, the gradient is large.

Image Divergent

In the sense of image processing, it is the inverse of gradient, this is a indicate of minor gray scale variation.In this aproach i used the channel Red and Green as vector field

Image Curl

Its a like a circular motion of variation of pixels. The curl operation measures the tendency of rotation around the point itself.


How can I manipulate the time signal to insert zeros in between frequency domain bins? I generate Multi tones in the frequency domain with constant amplitude and phase, for a fixed distance between tones.
As the figure below, where r is the change in frequency spacing between different symbols, n is the tone index starting on the left in the spectrum, N is the total amount of tones, M is the modulation order, and m is the symbol index.
fc=2.45GHz, f_delta= 1KHz, r=0.5,N=3,M=4.
I want to send the information between tones; I think I can do it by FFT frequency bins.
My question: How can I make some bins refer to one and other bins to zero between those tones?

Assuming the OP wants to manipulate a time domain function such that zeros are inserted in the frequency domain result, here is a simple approach:

Simply replicate the time domain samples and then divide by the number of repetitions to normalize and this will result in the same frequency domain result as the series that was replicated, with additional zeros inserted based on how many times the time domain sample was replicated.

Here is a simple example:

$fft([1, 2, 4, 2]) = [9,  -3,  1, 3]$

To insert one zero in between each sample repeat the time domain sample once and divide by 2:

$fft([1,2,4,2,1,2,4,2])/2 = [9, 0, -3,0, 1 ,0, 3]$

To insert two zeros in between each sample repeat the time domain sample twice and divide by 3:

$fft([1,2,4,2,1,2,4,2,1,2,4,2])/3 = [9, 0,0, -3,0,0, 1 ,0,0,3]$

It might be easier to understand how this occurs when you insert zeros in the time domain which causes repetition of the signal in the frequency domain over the duration of the DFT signal. This is explained in further detail here with regards to interpolation using zero-insert:

Upsampling and comb function


Uses of peaks in audio auto-correlation Do you know of any algorithms for audio processing that make use of peaks in auto-correlation function? I am currently involved in research on mapping methodology for dynamic algorithms. I came recently over an algorithm for detection of fundamental frequency in audio signals, which has a dynamic structure I am looking for - the execution time of the algorithm varies with the number of peaks in the signal's auto-correlation:

while (i < hsize)

    // Search for peaks in auto-correlation:
    if((r[i] > 0) && (r[i]>r[i-1]) && (r[i]>=r[i+1])){

            temp = 0.25f*(r[i+1]-r[i-1]);
            temp = r[i] - temp*temp/(0.5f*(r[i+1]+r[i-1]) - r[i]);
            n = n1[i] + n2[i];
            if ((rmax*n) < (temp*nmax)) {imax=i; rmax=temp; nmax=n;}.

            i++;}
    else{i++;}

This would be perfect as a studying example for my methodology, but the problem is that there is only a couple of lines inside the loop, so the overhead of adding my techniques becomes too big for the problem of this size. Now I am wondering if there exist any similar algorithms in audio processing, i.e. also making use of peaks in auto-correlation function but perform more extensive calculations for each peak?

The YIN Algorithm could be of interest to you. It's based on autocorrelation and is used for estimating the fundamental frequency.

A. de Cheveigné, H. Kawahara: YIN, a fundamental frequency estimator for speech and music, Proceedings of ICMC 2005. Download PDF.

A C++ implementation is also available.


Compressing plotted data with variable sample rate? I'm currently plotting high sample rate data by down-sampling by a factor and averaging. Most of the time this is working ok, but I am aware that at very low zoom levels (zoomed right out) the averaging is hiding vast amounts of data.

I have a timestamp for every data point, so are there any algorithms out there which will remove the pointless and 'less useful' points? I've done a bit of Googling and found papers retlated to time-series data compression, but that's not really what I'm after. I'm looking for something which will take advantage of the particular properties of plotting. ie. if you have 10 consecutive points with the same value, the intermediate ones can effectively be removed without losing any detail in the visual representation.

Here is what I did. The basic idea is to mimic what you would see if you were plotting every single sample, just do it much faster

Calculate roughly how many data points will fall in one horizontal pixel. There is no point in plotting more points than the are pixels on the screen.
Divvy up the data in blocks so that you roughly get one block per pixel
Calculate the min and the max for each block
For each block just add the max and the min in the correct order to your plot data

This way you get a plot array that has about twice the number of data points than there are pixels which is about the highest number that makes sense. You have to recalculate whenever you zoom, change the data range or re-size the window, but that tends to be fairly quick. If you don't want to handle window resizing, you can simply just calculate for the worst case (full screen).

Here is a code snippet that shows that in Matlab (it omits the ordering for brevity). The hi-res and the low-res plot look almost the same.

%% reduced sample plot example
n = 200000; % lots of samples
fs = 44100;
% noise with exponential decay
x = randn(n,1).*exp(-(0:n-1)'/fs/3);
% block approach
npix = 2000; % emulate an HD display
x1 = reshape(x,n/npix,npix); % chop into blocks
x2 = [min(x1); max(x1)];  % cal min and max
x2 = x2(:); % interleave
figure(1); plot(x); % plot with 200000 data points
figure(2); plot(x2); % plot with 2000 data points


Sampling Theorem: T-Sampled Signal I am stuck in a small question which I am trying to solve.

The signal $f$ defined by $f(t) := \text{sinc}(t)$ is to be sampled. How should I draw the illustrations of the $T$-sampled signal $f_{T} = (\mathbb{Z} \ni n\longmapsto f(nT)) $ for $T \in \left\{1,\frac{1}{2},2 \right\} $?

Thanks.

This is the same answer as I posted here:

Sampling Theorem

Sampling in the Time Domain created replications on the Fourier Domain.
The distance between the replications center is according to the sampling rate.

If you sample at $ T = \frac{1}{2 \Omega} $ than the replications are $ 2 \Omega $ apart.
The point here is to understand the replications are added to each other.
Hence if they are not far away enough they are summed (What we call aliasing).

Now just draw the Fourier Transform of you signal and according to the sampling rate add to it the replications.


Audio Signal Processing & DHT Is there any matlab (default) function for Discrete Hartley Transform? And is it applicable for audio signal? I want to apply it instead of DWT or FFT.

There is no in-build MATLAB function for DHT, but you can implement it very easily using the fft function:

function H=dht(x)
    X = fft(x);
    H = real(X) - imag(X);
end


Frequency response of a microphone using a sine sweep I want to determine the frequency response (magnitude, phase) of a microphone. I have another "good" reference microphone whose frequency response I know.

I understand that I can use a good speaker and measure the two mics simultaneously using white noise. Regarding the magnitude response, I can estimate the power-spectral densities by a method of averaged periodograms. Next, I can find the magnitude response of the device under test by subtracting the power spectral densities in logarithmic scaling and adding the known magnitude response of the reference microphone.

Averaging over periodograms seems necessary to control the variance of the periodogram estimator. But what about the phase response? Does a similar inconsistency exist for the angle of the fft also?

If I use an exponential sine sweep $x(t) \propto \sin(\omega_0 e^{a t} t)$, how do I obtain the frequency response? Do I average over the entire measurement time or do I do no averaging at all?

Edit:
I did not explicitly state but implied that I perform the measurements in an anechoic room.

Say $\mathbf{h}_1$ and $\mathbf{h}_2$ are the impulse response coefficients of mic1 and mic2, respectively, where $\mathbf{h}_1$ is known. Let $\mathbf{x}$ represent the input signal. Then, we can write the output of mic1 as

\begin{align}
\mathbf{b} = \mathbf{X} \mathbf{h}_1
\end{align}
where $\mathbf{X}$ is a Toeplitz matrix. If you want to find the impulse response coefficients of mic2, then you need to solve the following optimization problem
\begin{align}
\arg \min_{\mathbf{h}_2} ||\mathbf{X} \mathbf{h}_2 - \mathbf{b}||^2_2
\end{align}

The solution to the above problem is a simple least-squares solution. In the above formulation, I tried to express the convolution operation as a matrix multiplication. Find the filter coefficient response and express the transfer function in terms of $H_2(z)$. Once you know the transfer function, it is easy to determine the magnitude and phase response.


Mixed signals separation I have an unknown number of senors that sends the same pluse with different delays, the receiver receives a mixed signals. Is there any method to separate the signals and determine the number of sensors ?
Ps : I'm using MATLAB.

The typical method is to use the cross correlation function to recognize the start of each pulse sequence. Typically pulse patterns are chosen for good autocorrelation properties for this very purpose (likely the start of the packet as a header if there is a random data payload also involved). "Good autocorrelation properties" refer to an autocorrelation that is very strong when the sequence is aligned to itself but very weak for any other delays between the sequence and itself (white noise has this property). If the sequence correlates to delayed versions of itself, it will be difficult to differentiate in time when the sequence was received or to differentiate from separate sensors if they all use the same sequence.


Fit piecewise linear data What is a robust way to fit piecewise linear but noisy data?

I'm measuring a  signal, which consists of several almost linear segments. I'd like to atomatically fit several lines to the data to detect the transitions.

The dataset consists of a few thousand points, with 1-10 segments and I know the number of segments.

This is an example of what I'd like to do automatically.

I tried two approaches, naively (using only 3 segments). Surely there would be fancier methods out there.

 RANSAC, supposed to be a robust fitting mechanism. It's easy to stop the algorithm after a number of segments. However it may be difficult to enforce continuity between segments--as seems required in your application-- at least with a simple implementation.

As a proof of concept, I created an image from the data points so that I could use the RANSAC engine available in $ImageLines$, the line detection function of Mathematica.

 Fit a piecewise linear model using a general purpose minimizer. It's easy to enforce segments continuity. Interestingly, testing for residuals and other properties may provide enough information to determine automatically the number of segments--I've not tried it though. That's how it looks in Mathematica:


Frequency response and sampling theorem for triangular function The triangular function is defined as follows:
$h_l(x) = \begin{cases}1-|x|,&|x|<1;\\0&\text{otherwise}.\end{cases}$
According to ccrma.stanford.edu:
"If the output of the interpolator is also sampled, this can be modeled by sampling the continuous-time interpolation result, thereby aliasing the $\text{sinc}^2$ frequency response. [...] The Fourier transform of $h_l(nT/L)$ is the same function aliased on a block of size $ f_s=L/T$ Hz."
What I am thinking: If the Fourier transform $\hat{f}$ has bounded support (bandlimited), then $f$ has unbounded support (infinite length). This is a consequence of the uncertainty principle. But the triangular function (i.e. $f$) is $0$ for $|x| > 1$, so $f$ has bounded support (and $\hat{f}$ cannot be bandlimited). The sampling theorem requires $\hat{f}$ to be supported on some interval $[-B, B]$ (bandlimited). Then we cannot apply the sampling theorem, because $\hat{f}$ is not bandlimited. So aliasing will occur because we cannot find a good sampling rate.
Is my interpretation correct? Aliasing will always occur for $h_l$ no matter whether we look at $\text{DTFT}(\text{sample}(h_l))$ or at $\hat{h_l}$ (Fourier transform)? Then why is linear interpolation so popular in image processing (bilinear interpolation) if it has such a bad frequency response?

You're right that a triangular function is not band-limited and that it can't be sampled without introducing aliasing. However, this doesn't mean that linear interpolation can't be useful. The usefulness of linear interpolation depends on the data to be interpolated and on the desired interpolation factor. If the data have a low pass character, i.e., if they are already sufficiently oversampled, then linear interpolation can be good enough if the interpolation factor is not too large.
The obvious advantage of (bi-)linear interpolation is its simplicity and computational efficiency. In image processing it performs better than the common nearest-neighbor interpolation. However, if there are no real-time requirements and if computational efficiency is no big issue, then the standard interpolation method in image processing is the more complex bicubic interpolation.


Recover frequencies from IQ samples I have a set of IQ data. From those data I'm trying to get the amplitudes and the frequencies of my signal as I then want to plot them vs. time.

I am able to obtain the amplitudes by squaring both my I and Q values, summing them up and taking the square root of the sum; however, I am struggling to obtain the frequency.

I understand that I need to take a FFT to get from the time domain to the frequency domain but I'm not sure on what values I should apply the FFT. Should I do it individually on Is and j*Qs and then sum them up. Should I do it on the sum Is + j*Qs (which just gave me a new array of complex numbers of the form x + j*y). What role my center frequency / sampling frequency are going to play into that?

For context: I'm doing this using Python. (And I'm obviously pretty new to all of this.)

You should take the FFT on Is+j*Qs. If you have say 1000 IQ samples, you can form 1000 I+jQ complex values and take 1000 point FFT to get another set of 1000 complex numbers. If you plot magnitude of these output values, it will show you the frequency content of those 1000 IQ samples.

Since you took 1000 point FFT, each of those complex output values correspond to frequencies at multiples of $\frac{F_s}{1000}\text{Hz}$ where $F_s$ is the Sampling Frequency. For example, the first value will be corresponding to frequency at $0\text{Hz}$, next value is the frequency content at $\frac{F_s}{1000}\text{Hz}$, 3rd value corresponds to frequency $\frac{2F_s}{1000}\text{Hz}$ and so on..

If you are using numpy fft, the documentation should be a good enough guide to you.


Window function for signal with 3-5 sine waves inside Should I apply some window before FFT if I want to find 5 frequencies in signal? Signal is pure 5 sinusoids + some noise generated during recording on analog cable.

Windowing, in the context of DFT, is used to reduce the side-lobes or spectral leakage which happens because you work with finite length discrete signals.

Now, if is sum of sinusoids that are periodic in the observation time (that means that each sinusoid completes an integer amount of periods), then the DFT will be a series of delta function with no side-lobes.
However, in a more realistic scenario, the DFT of your signal should be effected by the spectral leakage. In that case, it is possible that some of the side-lobes will have greater magnitude than some of the sinusoids and if you try to select the 5 frequencies with the greatest magnitude, you will get the wrong results.
To sum it up simply: Yes, you should use some kind of windowing. It doesn't actually meter which window function you use. As suggested, Hamming window is a good choice.


Is it LTI system? given $H(z)$ is LTI Is it true or fale:

Given LTI, $H(z)$, with impulse response $h[n]$ that is not zero for all $n$. The following system is not LTI:

$$x[n]\longrightarrow\boxed{\downarrow2}\longrightarrow H(z)\longrightarrow\boxed{\uparrow 2}\longrightarrow y[n]$$

TRUE

Why? If I take delta to be the input, and then take another input: delta shifted by one...how I see that it will support the statement?

Lets say that $x_1[n]=\delta[n]$ and $x_2[n] = \delta[n-1]$, where $\delta[n]$ is Kronecker delta and that $y_1[n],y_2[n]$ are the output signals for $x_1[n],x_2[n]$.
The downsampling operator discards every other sample. Without the loss of generality, lets assume that we are keeping the even samples and throwing the odd ones. Therefore, after downsampling $x_2[n]$, it is all zeros, which means that $y_2[n] = 0$ for each $n$. Now since you say that $h[n]$ is not all zeros, it is guarantied that $y_1[n]$ will have some non-zeros values.
In other words, $y_1[n]\ne y_2[n]$ for some $n$, and you system is not TI.


Fragmented recovered image using a Phase Retrieval Algorithm named Error-Reduction After reading the paper Phase retrieval algorithms: a comparison by J.R.Fienup. I implemented the Error-Reduction Algorithm in MATLAB. The initial phase input was randomly generated, so the final recovered image in each execution is different. The only information given is the Fourier Modulus of the original image.
Sometimes I can get a pretty well result, but most of the time the result is a fragmented version of the original image as the figures below shown:

As you can see, the image of $\pi$ is "almost" recovered except that the head and the body of $\pi$ had been cut and swapped.
My questions are:
1. Is this a common problem using phase retrieval algorithm? What causes this?
2. Is there any method to avoid this?
3. Is there any other great methods in phase retrieval?

I think the effect you're seeing can be explained in terms of the shift theorem of Fourier transforms and the periodicity of the discrete Fourier transform.

It looks like the recovered object has been shifted upwards, and because the DFT is periodic, the top part of the object that's left the field of view has wrapped around to re-appear at the bottom of the field.
The reason the object is shifted is hard to see without knowing more about how your implementation of phase retrieval.  I can say that, due to the shift theorem, we know a shift in the object domain causes a linear phase in the Fourier domain, but does not affect the Fourier modulus.  So, the Fourier modulus alone can't help you find the correct shift.
I suspect you're leaving out some key information when you say "The only information given is the Fourier Modulus", because if that were strictly true you wouldn't be able to do phase retrieval at all.  As discussed in the paper, you need to have some other piece of information, typically the modulus of the object or the object support (the object support is a closed boundary that you know contains the object).  If you know the modulus of the object or at least its support, then reinforcing that constraint on every iteration of the error-reduction algorithm would prevent the object estimate from shifting.

As to your third question, I would point out that error-reduction is only the simplest algorithm discussed in the paper.  Gradient searches or the input-output methods described in the later sections are usually faster and/or more reliable.  In more recent years, successful methods that leverage nonlinear optimization algorithms such as BFGS have also become common. If you're interested, you might check out Prof. Fienup's list of publications.


Is it true that each and every signal has multiple frequency components? Usually when we observe spectrum of a signal,we see multiple frequency components, each having different contribution to amplitude

But is it true that each and every signal has multiple frequency components? For example common AC power has frequency 60 Hz, if we observe its spectrum, will we observe only one frequency components or multiple components?  Why if we observe multiple frequency components,despite the fact as we know beforehand that ac power has frequency 60 hz?

The only signal, that really has just one frequency component, is an infinite sine signal. Limiting the signal duration in any way is bound to produce other frequency components, as time limiting can be thought of as multiplying with a rectangle window, which translates to convolution with an si-function in the frequency domain, thus introducing new frequency components, however infinitesimal.

As for your example: the 60Hz AC power "signal" is messed up with all kinds of stuff coming from devices connected to the grid, that "backfire" into it. Also, the 60Hz are not really constant, they vary ever so slightly around 60Hz.


How can I decompose a signal into square waves? I am dealing with signals that are a superposition of different square waves with different amplitudes and phases. Normally, one would decompose a signal into sine waves with help of the Fourier transform, but in this particular case a decomposition into square waves would be much more effective. A Fourier transform would produce a very complicated spectrum, while a square wave decomposition should give just a few clear lines.

I know that such a decomposition is possible. In fact, I could use any periodic function as a basis for the decomposition and this is mentioned in many texts on the subject. But I could never find a formula or an explicit example for a decomposition into a non-sinusoidal basis.

My approach to decompose a signal consisting of the $N$ samples $x_k$, was to use a DFT-like formula
$$ u_k = \sum_{n=0}^{N-1} x_n \, \mathcal{R}_k(n)$$
where $\mathcal{R}_k$ is a real-valued square wave with a frequency $k$ times the base frequency. But this is certainly not complete, since I do not obtain any phase information for the constituent square waves, and I couldn't invert the procedure.

How can I decompose my signals into square waves with well-defined amplitude and phase?

What is described in the question is very near the Discrete Wavelet Transform (DWT) with the use of the Haar Wavelet.

The DWT decomposes a signal into a sum of dilated and translated orthogonal functions that do not necessarily have to be trigonometric. The DWT does not transform a signal from the time domain to a frequency domain but to a scale space where the "time" dimension is preserved. The Haar wavelet is effectively just one period of a square wave and due to its dilation and replication as the transform progresses it would appear as occurring at different frequencies. For more information on the link between decomposition level and frequency please see this link

Another transform that might be of help here, is the Walsh-Hadamard transform which does exactly that, decompose a signal into a sum of square waveforms which are orthogonal (Please note the sequency there as well).

For a brief example that seems to be near what you are after, please see this link

Hope this helps.


Thermography time evolution I have an imaging problem.  I essentially have a metallic plate with a number of holes through it.  I am looking at the number of holes with a thermal camera.  While looking at the plate I have two phases of this analysis, 1.a source of hot fluid is flowing through the holes and behind the plate 2.a source of cold fluid flowing behind the plate and through the holes.

During each phase I take a number of images N.  So in effect I have a set of random variables for each pixel in the image.  $H[x,y] = (x_1, x_2 .... x_n)$ and similarly $C[x,y]=(x_1, x_2 .... x_n)$ where H is the hot frames, C is the cold frames , x-y are the focal plane position of that particular pixel and 1,2,3 are the discrete time of the reading x.

As it stands I'm just taking the mean value for H and C then subtracting them to generate an image ($I=\bar{H} - \bar{C}$ where the bar denotes the expected value) to gain an image that has significantly decreased noise and removes much of the background effects.

This is really just to decrease the noise but I think there is more information there that I'm not utilizing considering this is not a wide sense stationary process, but is evolving with time as the temperatures reach equilibrium for each phase.

I'm looking for a way to utilize this data and for any recommendations on processing this time evolution to form an image.  Any suggestions?

I have tried Finite Difference Operator to approximate first and second derivatives, but this and various other algorithms have failed me.

I've had this post up with no takers so I will just put what I have found so far.

The algorithm that has resulted in the least noise for my application that I have tried is:

Let $X =[x1,x2,....,xn]$ be a sample of a function over discrete time intervals where x is a real number and represents the output of a particular pixel at $\delta t*i$ where $0<i<n$ .

Find a linear function $f$ that best approximates the data in the form (I tried polynomials of varying degrees but a linear function performed the best):

$f(x) = m *x + b$

Evaluate the function at ${t }= {{n*\delta t} \over{2}}$

So $f({{n*\delta t} \over{2}}) \approx Median $

My motivation was to find the median but it may be a somewhat bold assertion.  I wanted the median value of each pixel but the chances of picking a pixel due to noise was too great.  Plus because of the physical constraints I know that the data is either decreasing or increasing over the intervals with respect to the phase of data collection.

I do this for both hot and cold images to form the image using the notation above $I=\bar{H} - \bar{C}$.


How to estimate the frequency characteristics of a nonlinear chirp? Given a partial time-domain recording of a noise-corrupted chirp that asymptotically approaches a frequency, and the function that determines the non-linearity, how can I accurately estimate the final frequency it approaches?

Here's an example of what I'm talking about. A radio receiver hears noise until around sample x=5100, then the chirp starts at f=0, and quickly sweeps up to an unknown frequency, but then fades out around x=5500 to some annoying filter ringing caused by a hardware front-end band-pass filter.

X: sample count, Y: ADC count

Ideal solution will have the following characteristics:

Accurately predict the final frequency
Work over a range of final frequencies
Be invariant to different starting frequencies
Be invariant to different stopping frequencies
Relatively quick to compute

Approaches Tried So Far:

Estimating instantaneous frequency using zero-crossing, then curve-fitting. It seems that noise and filter-ringing causes errors in the frequency calculation that make it very sensitive to slight changes in waveform shape. It works better but not well if the signal is resampled to a higher sample-rate first.
Estimating instantaneous frequency using autocorrelation, then curve-fitting using the strongest lags from autocorrelation Function. This works well given knowledge of the incoming frequency so that the autocorrelation window size can be picked ahead of time, but fails when the ACF window isn't wide enough to see the whole signal. If the ACF is too wide, it'll sometimes see a harmonic as the strongest lag.
Curve-fitting the raw signal. This tends not to converge if there's even a little error in the initial conditions, or noise in the signal.
FFT/STFT, then curve-fitting the bins. What I ran into here was that it seemed the bins were destroying the shape of the signal. I couldn't find a good size of bin that gave enough frequency-resolution and time-resolution.

An example recording is available here.

How about just using $\dfrac{d\phi}{dt}$ to compute the instantaneous frequency?

GNU Octave code (I just made up a sample rate of 2 Msps):

pkg load signal;
x = csvread('nonlinear_freq_fit.csv');
y = hilbert(x);
z = y - mean(real(y));
dphi_dt = diff(unwrap(arg(z)));
Fs = 2e6;
f_inst = (Fs/2)/pi * dphi_dt;
k = [0:length(f_inst)-1];
t = k/Fs;
plot(k, f_inst);
axis([5000 5800 -2e5 2e5]);
title('Instantaneous Frequency');
xlabel('Sample index');
ylabel('Frequency (Hz)');
grid on;


Why $y[n] = x[-n]$ is not time-invariant? I followed these steps, but the answer still says that this system is time-invariant

let: $x_2[n] = x[n-k]$

$$\begin{align}
y_2[n] &= x_2[-n] \\
 &= x[-(n-k)] \\
 &= x[k-n] \\
\end{align}$$

and

$$\begin{align}
y[n-k] &= x[(-(n-k)] \\
 &= x[k-n]
\end{align}$$

where am I wrong in the analysis?

A time-invariant system is one that, when you shift the input signal, the output is shifted by the same amount.

A system that reverses the signal cannot be time-invariant because when you shift the input, the output is shifted the other way. $k$ and $-k$ are not the same amount.

$$
y[n-k] = x[k-n] = x[-n \mathbf{+} k]
$$


BIBO stable LTI system frequency response for this input signal? For an LTI system with bounded input and bounded output, I have the input $$x(t) = 5 + \cos(12t+\pi/4)$$ and output $$y(t) = 6\sin(12t)$$

It is said that the magnitude of the frequency response $|H(j12)| = 6$. I don't understand how this is true. Obviously the second term of $x(t)$ is scaled by $6$, but the first term needs to be scaled by $0$, right? Which means this LTI system isn't possible. Can anyone clarify?

For a sinusoidal input $x(t)=\sin(\omega_0t+\varphi)$, the output of a real-valued LTI system is given by

$$y(t)=|H(j\omega_0)|\sin(\omega_0t+\varphi+\arg\{H(j\omega_0)\})\tag{1}$$

where $H(j\omega)$ is the system's frequency response. The amplitude of the output is scaled by the magnitude of $H(j\omega)$ evaluated at $\omega=\omega_0$, and the phase is shifted by the argument of $H(j\omega_0)$.

If the input is a constant ($x(t)=c$), then in analogy with $(1)$, the output is given by

$$y(t)=|H(j0)|\cdot c\tag{2}$$

Since the system is linear, we can use superposition, and compute the output for each input component and then sum up all the individual components. So the output corresponding to the input $x(t)=c+\sin(\omega_0t+\varphi)$ is given by

$$y(t)=|H(j0)|\cdot c+|H(j\omega_0)|\sin(\omega_0t+\varphi+\arg\{H(j\omega_0)\})\tag{3}$$

So there is no contradiction between a gain of $6$ at frequency $\omega_0$ and a gain of $0$ at frequency $0$ (DC). From the given input and output you can deduce the values of $H(0)$ and $H(\omega_0)$.


Confusion about 'memoryless' meaning I am reading the book "Multiple input describng functions and nonlinear system design" written by A.Gelb and W.E.Vander Velde.

At some point it says:

  Single-valued characteristics are termed memoryless; multivalued characteristics are said to possess memory.

I am a bit confused about the meaning of the 'memoryless' term. I thought that memoryless meant a function that depends only on the value at that time instant and not at previous instants of time.

What am I missing?

The definition as you quoted it is only correct for static nonlinearities. The output of a static nonlinearity only depends on the input function directly, and not on its integral or derivatives. Otherwise, the nonlinearity is called dynamic (i.e. having memory). A static nonlinearity can also have memory if it is multi-valued, i.e. if there is more than one possible output value for a given input value (e.g. a system with hysteresis). Finally, a static nonlinearity is indeed memoryless if it is single-valued, i.e. if its output is uniquely defined by its input value:

$$y(t) = F(x(t))$$

where $y=F(x)$ is an injective function (i.e. each $x$ is mapped to only one value $y$).

So in sum you have

dynamic nonlinearity: has memory
static nonlinearity

multi-valued: has memory
single-valued: memoryless


Fir filter designing in matlab(1) Thanks porten!
Thanks Bitstream!

But our professor has said that design your own filter!(not used built in function of matlab).

I am designing my own filter.The problem is still there.First, i have to know filter's  order(N). One way is to find from its transition width. The pass-band edge frequency is 500Hz and stop band edge frequency is 700Hz. But i don't understand how to get it?
Can you help me?

There is no one perfect order "N".  You can always get a better (or at least as good) filter by increasing N, at the cost of more computational load.  So usually the question is "what is the lowest N that will be barely good enough?"

A relatively easy way to determine that is experimentally.  Guess what you think the optimal N will be, design a filter, see if it meets the filter criteria or not.  If it does with large margins, reduce N by a lot.  If it barely makes it either stop or reduce N by a little.  If it doesn't make the criteria, increase N.  Rinse and repeat.

frederic harris (for some reason he likes to have his name uncapitalized) gave the following formula in his book "Multirate Signal Processing for Communication Systems":

$ N \approx \frac{f_s}{\Delta f} \frac{Atten(dB) - 8}{14}$

where $N$ is the estimated filter order, $f_s$ is the sample frequency, $\Delta f$ is the transition band width, and $Atten(dB)$ is the amount of attenuation you want in the stop band, in dB.  I would not use the formula to get the optimal $N$, rather I would use it as a starting place for the iterative process described above.


frequency response of a signal in which two signal added in Here goes the problem I faced:

Sig1= A*sin(w*n)
Sig2= B*sin(w*n+phi)

Length of both signal is L.

Then, I put them in F(1:2:2*L)=Sig1 and F(2:2:2*L)=Sig2.

Now we know what would be the frequency response of Sig1 and Sig2 but is there any way we can find frequency response of F? FFT of sinusoid is the delta function.

Thanks

What you ask for should rather be called a frequency spectrum or Fourier transform of the resulting signal. The term frequency response is reserved for linear time invariant systems.

Anyway, the answer to your question can be found using the z-transform of the signals involved. First let's expand the two source signals into a series of unit delays:

$$S_1(z) = \sum_{k\in\mathbb{Z}} a_k z^k$$
$$S_2(z) = \sum_{k\in\mathbb{Z}} b_k z^k$$

Your combined signal can be written as

$$F(z) = \sum_{k\in\mathbb{Z}} a_k z^{2k} + \sum_{k\in\mathbb{Z}} b_k z^{2k+1}$$

and we can plug in $z=\exp(i\omega)$ to get the Fourier transform:

$$F(\omega) = \sum_{k\in\mathbb{Z}} a_k \exp(ki2\omega)+\sum_{k\in\mathbb{Z}} b_k \exp(ki2\omega)\exp(i\omega)$$

Using the Fourier transforms from the z-transforms above we can see that this actually equals:

$$F(\omega) = S_1(2\omega) + S_2(2\omega)\exp(i\omega)$$

You can interpret this result as the additive combination of the two original signal spectra, each scaled in frequency direction by $1/2$ and one phase shifted proportionally to $\omega$.

Since your original signals are of finite length, your initial spectra are sinc-shaped, centered at the same frequency and with a constant phase shift. I'll leave working out the details to you.


Best averaging to convert 50Hz Chi data to one hour averages? I have a $50\text{Hz}$ data set of temperature variance destruction or $\chi$ in the ocean for 365 days and hence is too bulky to handle. I want to convert it to hourly averages.

Simple question: How do I decide the best averaging method given that I do not have any pre-held notion about the data i.e I do not know what to expect? How can I justify its use for the current data? How do I make sure I do not lose the seasonal, diurnal trends etc. when I average?

My further calculations heavily depend on this processing and also this processed signal will give me spectra of different physics in the ocean. So I want to tread carefully. Any literature, suggestion should do to start off.

PS: In few posts and forums, I have read people commenting on moving averages as a "lousy" way to do it.

Essentially, you want to take your original signal and downsample it. Since you need to take care of aliasing before downsampling, you'd need to perform low-pass filtering of the signal.

Your original sampling frequency is $F_S=50Hz$, you want to downsample it to $1/h$ (1 per hour), which is $F_S/(60*60*50)$. So, your anti-aliasing filter needs to have cut-off frequency of $0.5/h$. After filtering you can just take every (50*60*60)th sample of the filtered sequence to have hourly averaged data. Note that downsampling to hourly averages requires, that the actual information in the signal does not vary too quick over time.

Code for some rudimentary filter design is below. Since your cutoff-frequency is so low, you need a very long filter to perform a good filtering. This is in tradeoff against calculation time.

Fs = 50
T = 24*60*60  # data of one day

cutoff = 0.5/(60*60)  # cutoff frequency

# Create some simulated measurement data
t = np.arange(0, T, 1/Fs)
x = np.sin(2*np.pi*t*3/(60*60*24))  # frequency: 3 periods/24 hours
y = x + 0.1*np.random.randn(*x.shape)

b = signal.firwin(15*3611, cutoff, nyq=Fs/2)
a=1

plt.subplot(1,2,1)
plt.plot(b)
plt.subplot(122)
w, h = signal.freqz(b, a, 32*4096)
plt.plot(w, 20*np.log10(abs(h)), '-x')
plt.xlim((0, 3/(60*60)))
plt.axvline(cutoff)

r = np.convolve(y, b)
r = r[:len(y)]

plt.plot(t, y)
plt.plot(t, r)
r_down = r[::50*60*60]
t_down = t[::50*60*60]

plt.stem(t_down, r_down)

As you can see, I get accurate downsampled data.


Why does a strange frequency appear in this image? Here is an image of a filter bank spectrograph, where the main signal is at 2.5kHz and the high signal is at 21.5khz,

What is the high frequency signal? I doubt that birds can make frequencies above 18kHz, so it is probably an effect from the recording equipment? Is it an issue from generic FIR BPF at high frequencies?

If the sample rate of data fed to a DFT or FFT for spectrum analysis was around 24 KHz, then what you are seeing is just the complex conjugate negative frequency image of the low (2.5 kHz) frequency audio data.  Basically, one should ignore any spectral information in a spectrogram for any frequency around or above half the sample rate, as it is redundant, just a (log distorted?) mirror image of the magnitude spectrum below.

If you want valid spectrographs for higher frequencies, use a higher sample rate, both for acquisition and analysis, with a higher anti-aliasing low pass filter cut-off.

But make sure to low pass filter any data to below half the sample rate before acquiring or decimating for any spectrum analysis, or you will see aliasing artifacts.


General question about DFT I am learning Fourier analysis and without any teacher, just trying to read books on my own. I think I have made some decent progress but they are a couple of points which are still very unclear to me and that I can't find explained in any of the books I have.

One of the sources that I found in this document, which is great because it teaches DFT without really using complex numbers. I shall say that I understand complex numbers and I am aware of Euler's Formula.

So what they say in this PDF/document is that in fact, in the simple case you can create a DFT by just using $N/2$ coefficients. I thought the choice of $N/2$ was related to the Nyquist frequency. If the signal contains $N$ samples then the Nyquist sampling theorem says that the signal can't contain a wave whose frequency is higher than half of the sampling frequency (hence the $N/2$ harmonics in the DFT).

So to me, explained that way everything made a lot of sense and in the simplest case you just needed to do something like this:

\begin{align}
a[k] &= \sum_{x = 0}^{N-1} s[x] \cos\left(2 \pi k {1 \over N } x\right), \quad \text{ for } k = \left\{ 1, 2, ..., \frac N2\right\},\\
b[k] &= \sum_{x = 0}^{N-1} s[x] \sin\left(2 \pi k {1 \over N } x\right), \quad \text{ for } k = \left\{ 1, 2, ..., \frac N2\right\}.
\end{align}

So this seemed simple. Now it says that when $k = 0$ and when $k = N/2$ then we need to divide $a$ and $b$ by $N$ or multiply them by $2/N$ otherwise. I understand why when $k = 0$, because it's the DC offset, but didn't really understand why you had to do the same thing when $k = N/2$ until I read this post

QUESTION 1: It seems to indicate that when you use the exponential form of the DFT then when $k = N/2$ then you have $\exp(\pi)$ which is equal to 1. Then it seems that in that situation the coefficient $a[N/2]$ has a particular meaning by I don't know which one?

Now this is where I am lost. In the "complete" equation for the DFT you don't compute $N/2$ coefficients by $N$ coefficients. That means that has soon as $k > N/2$ then the frequency of the harmonics is greater than the Nyquist frequency. I have illustrated this with the following image:

We have $N=8$ samples, thus the fundamental frequency is $\frac 18$ and we have harmonics: $1\cdot \frac 18, \ 2\cdot \frac 18, \ 3\cdot \frac 18, \ 4\cdot \frac 18$. However as soon as go above that then the harmonics go beyond the Nyquist frequency.
QUESTION 2: why do we test the signal with harmonics whose frequency go beyond the Nyquist frequency?

Finally, and I think this is actually related to question 2, I keep reading about positive and negative frequencies, but I just can make sense of this at all? This is my question 3.
QUESTION 3: could you please briefly explained if that's possible why we do speak and need positive and negative frequencies? Why do we need to care for negative frequencies?

The need for all these coefficients may not be obvious if you are looking at only real number input to the DFT.

But a complex vector input of length N will have twice the degrees of freedom compared to a strictly real input of that length.  The DFT, as an invertible linear transform therefore must therefore retain that same number of degrees of freedom, which for complex input ends up in the relationship between the positive and "negative" frequencies, or upper half of the transform, to produce complex exponential frequency results, instead of just strictly real cosine and sine functions.

For strictly real valued input, you don't really need to "care" about the other half, as it is duplicated information in the form of a complex conjugate (where the imaginary components of positive and "negative" frequency bins all cancel out to zero).  Thus the "mirroring" reduces the number of degrees of freedom to match that of the strictly real-valued input.  But one might still use the other half for computational purposes, as it makes the matrix math more regular and symmetric, and thus optimizable in the form of some FFT algorithms.

Thus the negative frequency coefficients don't need to be computed because they are redundant (unless the input is complex and not strictly real valued). In that case, the frequency coefficients above N/2 don't need to be computed because they are identical to the negative coefficients that don't need to be computed


How to take Processing/Scalloping losses in account in SNR/SFDR/THD/SINAD calculations? I'm in the process of developing DSP library for our systems.
I have developed "Windowing" functions (Hanning, Hamming, Blackman-Harris, ...).
Each of these windowing function has his own "Processing Loss" factor.
For instance, in case of Hamming window, I have found that such Processing Loss factor is 1.34.
How should I take this factor into account when computing SNR, SINAD, SFDR, THD ?

Some of the window performance metrics associated with the DFT are covered in great detail in:

F. J. Harris, "On the use of windows for harmonic analysis with the discrete Fourier transform," in Proceedings of the IEEE, vol. 66, no. 1, pp. 51-83, Jan. 1978.
doi: 10.1109/PROC.1978.10837
Abstract: This paper makes available a concise review of data windows and their affect on the detection of harmonic signals in the presence of broad-band noise, and in the presence of nearby strong harmonic interference. We also call attention to a number of common errors in the application of windows when used with the fast Fourier transform. This paper includes a comprehensive catalog of data windows along with their significant performance parameters from which the different windows can be compared. Finally, an example demonstrates the use and value of windows to resolve closely spaced harmonic signals characterized by large differences in amplitude.
keywords: {Discrete Fourier transforms;Fourier transforms;Frequency;Harmonic analysis;Oceans;Parameter estimation;Signal processing;Signal resolution;Signal sampling;Smoothing methods},
URL: http://ieeexplore.ieee.org/stamp/stamp.jsp?tp=&arnumber=1455106&isnumber=31261  Non-paywall version here.

Not all the metrics you ask for are covered, but there are several which you haven't identified that are discussed. Perhaps if you said more about who might be interested in using your library, more specific suggestions are likely to be offered.


Computing the Power Spectrum Density (PSD) on a CSV File in Python im quit new with signal processing and im trying to calculate the PSD of a signal im sampling. the signal is an output of a DC buck converter

this is the code im using and this is the plot im getting

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

filename = 'scope_yos2.csv'
data = pd.read_csv(filename)
ConvertToMatrix = data.values
time = np.delete(ConvertToMatrix,[1,2],axis=1)
voltage = np.delete(ConvertToMatrix,[0,2],axis=1)
NumOfSampels= len(voltage)
plt.plot(voltage)
voltage1 = voltage.transpose()
samplFreq = 1.25e9
Pxx, freqs = plt.psd(voltage1,
NFFT=256,Fs=samplFreq,
detrend=mlab.detrend_mean,
window=mlab.window_hanning,noverlap=0,sides='onesided',scale_by_freq=True,
return_line=None)
plt.show()

first image is the DCDC ouput
second one is the PSD plot

the csv file contains two columns first is the time ax and the second one is the amplitude

how come the plot is so smooth. this is the plot im getting using the welch algoenter image description here

also theres a argument called **kwargs , what is used for ? , as you can see im not using it

thanks

I have not used the mlab functions in matplotlib, but from what I can discern from your time plot, I see a pattern that repeats roughly every 15,000 samples, at 1.25 Gsps according to your code.

This means we should see a strong spectral component at $1.25e9/15000=83.3$ KHz.

The reason we do not see this, I suspect is the relatively short FFT length (or very high sampling rate compared to the spectral occupancy, or both!).

Let me explain: The FFT length is currently 256 samples. At 1.25GHz this is only $256/1.25e9 = 204.8$ ns.

The frequency resolution prior to using any additional window (meaning using a rectangular window) is $1/T$ where T is the length of the FFT in time (for further detail on that see:
What happens when N increases in N-point DFT )

So in your case the frequency resolution is greater than $1/204.8e-9$ = 4.88 MHz!!  Any window will reduce sidelobes but at the expense of increasing the mainlobe, so the width of your mainlobe given the Hanning window used will cause the resolution to be even wider than 4.88 MHz. (But definitely good to include windowing, otherwise the additional sidelobes would be much higher). You can see this post for additional details on the windowing choices:Why would one use a Hann or Bartlett window?  (From this we see the Hanning window is approximately twice as wide as the rectangular window, so your example would have approximately 10 MHz of resolution bandwidth per bin in the FFT)

Looking at your plot, it looks like we see exactly this, which is the FFT bin of your strong DC component, and this single FFT bin is approximately 10MHz wide in bandwidth! So that is burying any signals of interest that you want to see.

My recommendation is to increase the FFT length within what your processing will allow. And if that starts to get problematic then reduce the sampling rate to get further resolution (however there are challenges in doing that properly that I won't get into except to say that you must be sure higher frequency components are filtered out prior to resampling). Without requiring a sampling rate change, you would want to use an FFT length of at least $1.25e9/10e3$ = 125,000 to get 20 KHz of resolution BW while using the Hanning window. That's a reasonable length and I suspect you can go much higher than that for finer resolution.

Also as a side note I don't see where the impedance is set, so I guess it is assumed to be normalized to 1 ohm. Just mentioning in case you cared about knowing the actual power levels.


Phase changes of BPSK modulated carrier after moving to baseband Given a BPSK-modulated carrier what would happen to the phase changes in the signal after the carrier is moved to baseband?

Would the time intervals between the phase transitions in the time domain remain the same?

There is no change - the complex envelope of the signal (magnitude and phase) at any given carrier frequency is no different from when that carrier frequency = 0 (DC, at baseband). To hold this equivalence the baseband signal must be complex (often represented with I and Q data paths) unless the signal is completely on the real axis... for example If the BPSK Signal at the carrier has a 10 degree phase offset (such that the carrier frequency is toggling between 10 degrees and 190 degrees), the baseband signal will do the same as long as it is allowed to be complex.

For this reason it usually doesn’t make sense to actually model a sinusoidal carrier when doing simulations but work directly with the complex baseband signal.

The further details given on upconversion and downconversion in this post may help provide further insight:
Frequency shifting of a quadrature mixed signal


Confusions in windows based design of FIR filters? I am reading the proakis book"DSP USING MATLAB 3rd Edition"
In chap 7, article 7.3, i came across example 7.8,where i saw formulas for pass band ripple and stop band attenuation
I have highlighted them in attached snapshot but i couldn't find out how they were calculated
Please also include reference of these formulas of Rp and As, into your answer, as i could not find their reference in book
and also what is "delta_w"?Is it "transition width"?

From the code you can see that they compute the maximum passband ripple by computing the minimum value of the filter's magnitude response in the passband. Similarly, they compute the maximum value of the magnitude response in the stopband, corresponding to the maximum stopband ripple.


Demodulation with MATLAB I am trying to listen VLF radio signals. I have a recorded wave file (download) and here is the frequency spectrum:

I have produced this spectrum using this matlab code:

signal= wavread('/Users/ecabuk/Downloads/DR0000_0165.wav');
X_mags = abs(fft(signal));
bin_vals = [0 : N-1];
fax_Hz = bin_vals*fs/N;
N_2 = ceil(N/2);
plot(fax_Hz(1:N_2), 10*log10(X_mags(1:N_2)))
xlabel('Frequency (Hz)')
ylabel('Magnitude (dB)');
title('Single-sided Magnitude spectrum (Hertz)');
axis tight

It has two clear picks near at ~24KHz and ~25KHz.

How can I crop 24KHz and then return to time domain and save as a wave file, so I can hear that. (and same thing for 25KHz of course)

I am not sure but am I trying to make a Demodulation process?

You are asking 2 questions:

How to extract those labeled frequencies and remove the rest of the data?
How to do Demodulation?

Well, regarding your first question, a dirty solution would be by working the the DFT data.
Just zero all frequencies but those you're after and apply Inverse DFT (ifft).
Leave some margin around your wanted frequencies.

About the second question, look here:
http://en.wikipedia.org/wiki/Demodulation

Basically, in order to shift data in the frequency domain you multiply by pure harmonic signal in the time domain.
You chose the frequency of this harmonic signal s.t. (At least this is the classic case) the frequency you want will move to the DC frequency.
Then you apply Low Pass Filter to leave only the data you wanted.

This is the opposite of the Modulation process.


Algorithm for identifying beat grid in music with non-constant tempo Are there any algorithms for automatically fitting a beat grid to a song (tempo matching) when the song has a slowly changing tempo or BPM?

The output would be a vector of onsets.

Here's a picture of fitting a grid to a waveform for a part of a song with a constant BPM:

I briefly worked on tempo and beat estimation for electronic dance music and used a couple of correlation techniques for this. Basically, I used a combination of auto correlation and cross correlation using a grid technique to find the tempo and then used the tempo to find the beats. This was for a constant tempo though. I'm not sure how well that would work on non-constant tempo. You can take a look at my code and detailed algorithm on my Github repo here.

As per the latest MIREX results (a worldwide competition for some music information retrieval tasks), you should check out the madmom library. They seem to have amazing results and use machine learning for their approaches.

What is your exact use case? Let me know if you find something useful.