Spread spectrum receiver I want to design a receiver for a spread spectrum system. I know that two Pseudo sequence with two different length must be used in the receiver, but I can' t manage this with common receivers(such as integrate and dump and so on...) , can this be done?
Regards.

Just despread the same way with two different integrate and dump times for your two symbol types with appropriate blanking on each between longer and shorter code sequence periods.

I would be tempted to do all the acquisition and tracking on the longer symbol and then simply despread the shorter symbols using the synchronization from the longer one— with that approach envision a receiver for the just the longer sequence with 62 chip gaps between each symbol.

Alternatively if each code repeats at the same location in each symbol, given the relatively short sequences, a sliding complex correlator using I and Q FIR filters with the sequence as the coefficients (in reverse order) using two samples per chip would be a very simple approach —- the complex I and Q output of the filter would have a spike in time after every alignment with those symbols (0 or 180 degrees once aligned depending on the data symbol and given I and Q would have phase to provide for carrier recovery).  Detect the spike with a threshold detector on the magnitude for acquisition, track the complex peak for carrier tracking and track the two samples on either side of the peak for early-late timing tracking. The peak sample will give you the same performance as an integrate and dump of the multiplication of the received sequence with the output of a code generator but given your peculiar waveform structure this may be a lot simpler.

Once all timing is established, a single 31 tap FIR filter on I can simply demodulate the other two symbols.

Summarizing this FIR approach:

One 126 tap FIR filter on I  for carrier and timing recovery/tracking and demod of 63 chip symbols.

One 126 tap FIR filter on Q for carrier and timing recovery/tracking.

One 31 tap FIR filter on I for demod of the 31 chip symbols.

As @WilliamHa points out in the comments, this approach will only work if the same code sequence appears in each symbol, starting at the same location each time (beyond being inverted to designate a "1" or "0" data bit), but will not work if the code generators run continuously causing the code starting location to shift by one sample at the start of each longer sequence, and at the start of the pair of shorter sequences, due to the relationship between the two codes.  In this case I would approach it using continuously running correlators where the outputs are blanked during the other symbols, using the longer code for tracking and carrier recovery and using the synchronized clock from that to demodulate the shorter codes. Similar to the FIR filter filter approach, there would be 6 correlators for the longer code: 2 for Early (I and Q), 2 for prompt (I and Q) and 2 for late (I and Q). Each of these should be running at at least 2 samples per chip. Only one I correlator is needed for the shorter codes at one sample per chip as the carrier and time reference is established from the longer code. The added complexity will be in timing the integrate and dump outputs of the correlator such that they are blanked during the shorter code sequences, but this would be determined in the timing recover process such that accumulation starts at the estimated start of the longer code and dumps after 63 chips, and then waits 62 chips before accumulating again. This structure assumes the code generators run continuously with the same chip rate, and the correlator structure is done by multiplying the code generator output with the received signal after having been down-converted to baseband. With this approach Carrier Recovery is done by monitoring the rotation of the I and Q prompt symbols of the longer code; a carrier error is determined using a complex conjugate multiply of one symbol with the previous and then multiplying that time I prompt to strip the +1/-1 data. Clock Recovery can be done by monitoring the magnitude of the Early and Late outputs of the longer code, where clock error = Early-Late. Acquistion is done by setting a threshold on the magnitude of the longer code.

Using the cross correlation property of the FFT, all the correlators could be done with FFT's as well, including one-shot joint frequency and delay acquisition if processing allows. I explain that approach further in this post: GPS signal acquisition


looking for simplified version of motion capture for simplified situation where objects are solid and have no joint (no internal degree of freedom) The whole question is about a simplified algorithm for tracking a single solid object which has no joint, with assumption that i have full control on adding some markers in real world to it.

For example i want to find full coordinates (Center of orange Coordinate system and direction of each orange axis) of the car in the following image:

Assume i have control over adding as many markers(Stars) with required features(Shape, color, position) to my object(dreamy car) in real world.

Notice i am looking for some algorithms which mainly rely on finding position of my markers, which could be detected very fast by some simple color threshold.

Fast way could be to convert your image to HSV format and then apply color thresholding as you have mentioned. Here is a short snippet which could get you started quickly.

import cv2 as cv2
import numpy as np

fn = 'image_or_videoframe'
# OpenCV reads image with BGR format
img = cv2.imread(fn)
# Convert to HSV format
img_hsv = cv2.cvtColor(img, cv2.COLOR_BGR2HSV)

# Choose the values based on the color on the point/mark
lower_red = np.array([0, 50, 50])
upper_red = np.array([10, 255, 255])
mask = cv2.inRange(img_hsv, lower_red, upper_red)

# Bitwise-AND mask and original image
masked_red = cv2.bitwise_and(img, img, mask=mask)

in this case, the red is filtered in the image and masked_red would contain only the red pixels in the image.

After the filtering, depending on the requirement you can use houghlines to detect the lines. It is not trivial but it is definitely doable.

Hope this helps and all the best :)


How to correct phase of Diagonal Volterra kernels obtained with exponential sweep sine? I'm trying to model a non-linear system using non-linear convolution with Novak's (2010) synchronized exponential sine weep (SESS) that models them with a Generalized Hammerstein (Volterra diagonal).

I get pretty clear kernels in terms of SNR but clearly phase-distorted. It's a common problem in real estimation.

I can theoretically correct them with a time compacting inverse Kirkeby filter (Farina, 2007).

[EDITED AFTER FIRST ANSWER]
I use Novak's himself code from http://ant-novak.com/Swept_Sine_Simulation.php

Changes: Doubled both number of kernels and sampling frequency and raised a lot SESS duration for better kernels separation.

fs = 192000;         % sampling frequency
f1 = 10;            % start frequency
f2 = 12000;         % final frequency
T_ = 50;            % approximative time duration
N = 8;              % order of nonlinearities

L = round(f1/log(f2/f1)*T_)/f1;
T = L*log(f2/f1); % this definitive sweep duration in theory makes Novak's SESS more robust to phase issues (synchronized)
t = (0:ceil(fs*T)-1)./fs;
x = sin(2*pi*f1*L*(exp(t/L)-1)); % sweep

Novak uses an elegant frequency domain definition alternative of the regular inverse sweep with amplitude correction. I've tested both with similar results, although Novak's frequency method gives better SNR.

The other change is that I model a "real" system, well, only a guitar amplifier VST plugin. So I use the plugin's output to the SESS for the deconvolution stage.

I'm not 100% sure, but the kernels look extremely "right-sided" like in Farina's 2007 paper.
The frequency magnitude of the kernels is coherent with a high gain guitar amplifier VST (even harmonics stronger). The kernels phase is a mess, triangle-saw like.
When I test the model with a dry guitar sound as input the output barely improves a regular (first kernel) equalization matching.

Both in the SESS (plugin input/output) and simulation (dry guitar signal) stages I apply a steep bandpass in the range 100hz-10Khz. The plugin is configured to do x2 oversampling (384Khz) to avoid aliasing artifacts.

A potential pitfall in non-linear modeling is using different signal levels. The sweep test needs to be carried out at a sufficiently high signal level to cover the maximum possible real-world input levels, but not much higher as this will require more kernels for the same effective quality. Usually real-world signals to the system are much of the time at lower level than their peak level, with the 1st order linear kernel dominating.

To me your kernels look just fine. I don't know how common it is to see that a higher-order kernel (3rd order in your plots) is more prominent than the 1st order kernel. The saw-like phase comes from the large initial delay, which gives a steep linear dependency of phase on frequency. Wrapped phase of $\pi$ is the same as a wrapped phase of $-\pi$, and you are seeing wrapped phase, which gives the sawtooth pattern. If you remove the initial delay in the inspection of the results, you should be able to unwrap the phase more easily for inspection.

Most systems you would model by the technique will have right-sided kernels. The systems are causal and cannot store energy for very long after they have received an impulse as input. That would require a delay line. Instead, much of the energy passes through almost instantly as the initial spike, and the residual energy is released as the tail.

I suggest you test with a system you have complete control over, for example a memoryless non-linearity that is the sum of the $N$ first Chebyshev polynomials of the first kind, one for each of the $N$ harmonics you want to generate, including the 0th and 1st harmonics in the count. Then any phase shifts due to errors in the algorithms/software should be clearly visible. If the algorithms work correctly, each harmonic impulse response should be a delayed sinc, and the the coefficients of the memoryless nonlinearity polynomial should be resolved as the amplitudes of the scaled impulse-like kernels.

What I have done in my own implementation is to generate the sweep and its fade in and fade out in the frequency domain by programming the timing into the group delay, and amplitude into the magnitudes, rather than using a direct time domain equation for the sweep. Either way should work and will have its own problems. The frequency-domain approach gives amplitude oscillations near the fades. In the frequency domain it is easy to correct the non-whiteness of the sweep in the inverse filter. (I can't publish my implementation in detail at the moment.)

If there are phase problems, I'd try $\cos$ instead of $\sin$ in the time-domain sweep formula, and fade the sweep in and out. The fades will take care of eliminating transients at the ends. The initial and final phases of the sweeped sinusoid could be made such that the value of the sweep is zero both in the beginning and in the end, but with the fades in and out that is no longer needed for transient reduction, because the fades gently bring the value to zero in any case.

I don't see a problem with the sweep going all the way to the Nyquist frequency, if the non-linearity is analog and the analog-to-digital converter (ADC) has an anti-alias low-pass filter at the Nyquist frequency. This way the harmonic impulse responses and consequently the Hammerstein kernels have only non-aliased frequencies up to Nyquist frequency. It is sufficient information to model the non-linear system in the sub- Nyquist frequency band if the system is compatible with the Hammerstein model. Only the implementation has to be careful not to introduce aliasing. In any case, you'd typically want the linear kernel to have frequency response information up to or at least near the Nyquist frequency, which dictates the sweep end frequency.


How does kernel work in unsharp masking? I applied the kernel below to the picture for unsharp masking

k = [-1 -1 -1; -1 8 -1; -1 -1 -1]

I got !

Then I changed 8 to 14. And I got
How does kernel work for unsharp masking the original picture and what does exactly do the number in the middle ?

The first kernel that you gave,

$$
\left[ \begin{array}{ccc} -1 &&-1 &&-1 \\
-1 && 8 && -1 \\
-1 &&-1 &&-1 \end{array} \right]
$$

is interpreted as follows: for each pixel in the image, take 8 times its value (hence the 8 in the center of the kernel) and then subtract the values of all the 8 neighboring pixels (hence all of the surrounding -1 values). This is a form of a highpass filter, which tends to accentuate edges.

Think about what you get if you were to run this kernel over a patch of a constant color. If the region has a uniform color, then it doesn't contain any edges. The kernel above would yield a value of zero when applied to a 3-by-3 pixel block that has the same color, which is desirable. Not only is the filter highpass, it also completely attenuates zero frequency, which is also good.

If you change the center value to 14 instead of 8, it no longer has the property where it will output zero when applied to a region of constant color. Essentially, the kernel will have a different spatial frequency response, which leads to the differences in the second output image.


Confused about HSI and HSL color spaces I am reading a book for image processing algorithms (Andreas Koschan, Mongi Abidi: Digital Color Image Processing) and for the contrast algorithm it says that I can either go RGB->HSL or RGB->HSI first and than apply a contrast technique for grayscale images, on the lightness component.

Then it gives this formula only, not other formulas for color conversion:

$$L(x,y) = 0.299*R(x,y) + 0.587*G(x,y) + 0.114*B(x,y)$$

This formula is neither for the L in HSL, neither for the I in HSI and that is what confuses me.

Can you explain what this formula stands for, and maybe help with the appropriate conversion formulas?

Thanks

L(x,y) here is a function for luma, a non-linear approximation of luminance. As you can see, it's the weighted sum of the R'G'B' values. The prime symbol (') indicates that they have been gamma-corrected for perceptual uniformity. The coefficients used in this case are standardized by ITU-R Recommendation BT.601. Other standards, such as BT.709 (which may be preferred for its higher accuracy), exist as well.

See: http://en.wikipedia.org/wiki/HSL_and_HSV#Lightness. I can't answer whether you can use luma and still say that it's a component of HSL or HSI. If the author of the book gave you that formula, I suspect you can assume so, at least for that text.


Audio output quality When a digital device (PC, portable media player, etc.) plays a digital audio file (ogg, mp3, flac, etc.) the audio signal output is always the same, regardless of the device type/brand, am I right? Thus the sound quality should not differ. For instance, I have an iPod that plays an mp3 file. If I play the same file on a different brand media player the sound quality should be the same since it is the very same signal (no equalizers, or any sound alterations).

Is this true? If it is true then I reckon it is only the headphones/speakers that matter as far as the sound quality is concerned.

P.S: all of the answers are really nice! wish i could accept 'em all

Unfortunately, the real-world gets in the way - there's a whole host of things between you and the "bit-perfect" waveform that comes from the MP3 decoder:

As you noted, user controllable filtering (graphic eq for example) - this may be disabled.
"improvements" to the sound made in the digital domain (for example compression to make it sound "louder" and bass boosting) - this may also be disabled
The internal processing is often done at a greater bit depth than the output DAC capability - for example, 32 bit processing but only a 16 bit DAC.  There must then be a process to reduce the bit depth.  This may or may not include dithering of some sort.  The results of this will sound different depending on the implementation.
The digital to analogue converters may be of differing quality, which will impact on (amongst other things) the linearity and noise floor of the analogue output.
The clock feeding the DAC will also have an impact on the performance - higher jitter will be noisier.
The headphone driver will also have various analogue characteristics - eg. noise floor, frequency colouration etc.


what Graphical DSP programming tools are there for ASIO4ALL? Many years ago I had a PCI sound card (I think it was some kind of SoundBlaster but not sure, this was maybe 15 years ago) with a 5 1/4 disk sized front panel with various ins and out.
I remember using ASIO4All but there was also a graphical component where I could program the DSP directly by creating a visual circuit and add effects and other features, and I used this to create a multitrack studio with talkback capabilities to the artists.
But I can't for the life of me remember what it was, or find anything similar on the webs.
is the 10k DSP chip still a thing? If not, what is the modern-day equivalent, and is the visual programming still available for it.
Thanks in Advance.

ASIO4ALL is a hardware independent low latency ASIO driver for WDM audio devices. As being a (wrapper) driver it does not support programming other but change of settings (device/channels selection, bit depth, sample rate, I/O buffers (latency)) through control software.
Regarding DSP programming tools, could you mean Creative/SoundBlaster specific kxproject driver/software?
Edit 1:
Creative/SoundBlaster DSP chips are listed here.
There are also software based methods as well as like EqualizerAPO  (system-wide equalizer for Windows) which lets you route channels, EQ filtering, convolution etc..


Manual Extraction of Filter Phase Response with MATLAB using sweep sine I'm currently at a loss for why my code seems to not be working.  I'm doing a sine sweep on a filter to calculate its freq and phase response.  The frequency response checks out, however the phase response doesn't look right when compared to MATLABs freqz function.  Plus since I know the filter is a comb, I can tell its wrong as well.

Here's what I've got so far.

%   Filter Coefficients
b = [1];
a = [1, 0, 0, 0, 0, 0, 0, 0, 0, -0.96];

%  sweep is my sine wave that goes from 0->fs/2
filtered_sweep = filter(b, a, sweep);

% This section correctly Grabs and plots the Freq Response
fft_sweep = fft(filtered_sweep);
freq_response = abs(fft_sweep);
freq_response = 20*log10(freq_response/1);      % Convert to dB scale
freq_response = freq_response(1:end/2);         % Cut to pi fs/2
figure;
plot(freq_response);

And Now here is where I try and get the Phase two different ways with the same wrong result.  From what I understand to get the phase you just need the angle of the complex values.  In MATLAB I thought angle() should work, but just to be safe I tried to manually compute the angle myself using atan with the imag and real components.

% Pre Allocate Vector
% Manual Angle Extraction without using Angle()
phase = zeros(1, length(fft_sweep));
for i=1:length(fft_sweep)
    phase(i) = atan2(imag(fft_sweep(i)), real(fft_sweep(i)));
end

% Using angle()
phase_response = angle(fft_sweep);
phase_response = phase_response(1:end/2);
figure;
plot(unwrap(phase_response));
grid on;

Both methods give the same incorrect phase response.  Any help would be appreciated.

You missing one, very important step. In order to get the impulse response of your filter while using the sweep sine, you must get the output of your system $h[n]$ to excitation with sweep-sine signal $x[n]$:

$y[n]=h[n]\star x[n]$

In your case that's the filtering with filter function. Now in order to get the impulse response you must convolve your signal with the inverse filter $f[n]$. Generally that's the time-inverted sweep signal (although special kind of pre-conditioning, amplitude scaling, etc. might be done as well):

$h[n]=y[n]\star f[n]$

Now you have your impulse response that you can use for further analysis. Please take a look at this answer tackling very similar problem.

For more theory please see:

  A. Farina - Advancements in Impulse Response Measurements by Sine
  Sweeps

  A. Farina - Simultaneous Measurement of Impulse Response and Distortion with a
  Swept-Sine
  Technique


phase difference detection For application of syncing sampling clock phase to carrier phase of received signal where the received signal is sampled with 10 times the carrier phase and modulated with the PSK shown in the figure (Baud rate up to 0.5*fc)

I had many ideas to apply on my system (where interpolation is not possible) and I need your help figuring out which is the best:

Timing Error Delay: most proposed algorithms require interpolation
Hilbert filter:where generating I/Q requires shifting 2.5 samples or Fractional Delay Filter (where I need help with determining the best FIR implementation of such a filter)
thanks in advance

For fractional delay filters see this article: https://www.dsprelated.com/showarticle/22.php.

Timing error delay is not used for carrier recovery. Carrier offset and symbol timing offset are not the same thing.

Carrier recovery can be done at 1 sample per symbol assuming you do timing recovery first. For this a simple phase detector is the cross product phase detector where you derive phase from the imaginary term of the complex conjugate product.

The Gardner Timing Error Detector is a good choice as it can determine timing offset over a relatively large carrier offset. If resolving timing offsets to within 1/10th of a symbol then this would be a very simple approach without further interpolation. However, I question why the OP precludes interpolation; combining a timing error detector with a polyphase interpolator is an excellent approach to resolve timing with fractional delay filters without having to increase the sampling rate- however interpolation techniques are still applied.

I have more information in other posts for both of these approaches listed below:

Gardner TED:

Isn't Gardner's algorithm and Early-Late gate the same thing?

Fractional spaced equalizer + timing (clock) recovery

Phase Detectors and Carrier Recovery Loops

Phase synchronization in BPSK

More Carrier Recovery and Timing Recovery

How to demodulation BPSK in real case(through a channel)?


Minimum recording time for a signal I am familiar with both Nyquist frequency and Nyquist rate. What I cannot seem to find information on is the minimum time necessary to detect a signal.

In my use case I am using a multiplexer to measure signals from multiple analog sources, each source carrying information from 60 to 100 Hz. I would, ideally, like to take these samples as quickly as possible.

Lets assume that my lowest frequency is 10 Hz. to detect a 10 Hz signal, you need at least 200ms of data, right? You cant detect a signal that comes once every 100 ms without measuring for at least twice as long. Or can you? This is part one of my question.

Part two is, how does the signal appear if I tried to measure a 10 Hz signal using only a 100 ms bin? Obviously it appears as a random point, but does it mask itself as a higher frequency? Is it just combobulated noise? What does the aliasing look like?

Lets assume that my lowest frequency is 10 Hz. to detect a 10 Hz signal, you need at least 200ms of data, right? You cant detect a signal that comes once every 100 ms without measuring for at least twice as long. Or can you? This is part one of my question.

Consider the Discrete Fourier Transform (DFT) of that "subsampled" block:

$$ X_k = \frac{1}{N} \cdot \sum_{n=0}^{N-1}x_n \cdot \left( \cos(\frac{-2\pi kn}{N}) + i \sin(\frac{-2 \pi kn}{N})  \right)$$

This formula has "no problem" in producing the $k^{th}$ sum, when $k \ll N$. So, yes, information about that harmonic is still available in the signal.

At the same time, notice that now, this sum is going to be scaled by the limited availability of data. So, if you "subsample" a 10Hz sinusoid by examining a block of data that contains half a period, you can still evaluate the sum for that low frequency, but the amplitude that you will "read" back will be half the true amplitude, simply because the sum did not have enough samples to "evolve" in.

In a noise free environment, that would not be too much of a problem but in practical conditions, this extra scaling would make reliable detection / extraction of that harmonic even more difficult in a noisy environment.

  ...how does the signal appear if I tried to measure a 10 Hz signal using only a 100 ms bin? Obviously it appears as a random point, but does it mask itself as a higher frequency? Is it just combobulated noise? What does the aliasing look like?

I am not sure how does the "...random point..." represent the signal in here (?), but the closest you can get to how that would look like is to take a "window" of the samples of a 10Hz sinusoid. Obviously, if you are sampling this at different (or random) times, you will be "catching" it at some random phase. So, the amplitude component will be more stable than its phase.

This "mental picture" probably takes us nicely into the discussion about aliasing:

Strictly speaking, there is no aliasing here because $F_s>f_{foi}$ (frequency of interest).

BUT!, if you are trying to estimate the amplitude and phase of those components via this kind of "blocked" recordings, then the frequency of sampling of the block does come into play.

It is not exactly clear form your question if the frequency by which you collect blocks of data is stable or even if, through the analog muxer, the signals are truly sampled at the same time.

So, to cut a long story short: If you are sampling in windows but do it regularly, there should be no aliasing. Your estimation of the amplitude will be scaled and the estimation of the phase will have a systematic error that is proportional to the rate at which the blocks are acquired and the initial phase that the first block happened to "catch" the sinusoid at the input. But, if you sample these blocks irregularly, with some jitter because of whatever reason, then this will creep in to your estimates. Again, the amplitude will probably be more stable than the phase but the latter will now vary in an irregular way that depends on the jitter.

Hope this helps.


Fourier Transform of Morlet wavelet Function? One definition of the Morlet wavelet function is given by:
$$\frac{1}{\sqrt{\pi f_b}}e^{\frac{-t^2}{f_b}}e^{j2\pi f_c t}$$
The Fourier transform of this equation is:  $e^{-\pi^2 f_b(f-f_c)^2}$ (is it right)?
First I attempted to plot the FFT of Morlet function by FFT function in Matlab,  then I've plotted the Fourier transform function directly. I've expected to see same plots, but unfortunately they have different magnitude (why?). I publish the plot and my codes, do you have any idea why this happened? Is there any problem in my codes? Or any problem in the obtained Fourier transform? Or any problem in Matlab built-in FFT function?

I wrote this script to get those plots:
clear all;
fS = 500;
tStart= -4;
tStop= 4;
timeVector = linspace(tStart,tStop, (tStop-tStart)*fS );
fC = 2;
fB=2;
timeMask = zeros(1,length(timeVector));
timeMask((timeVector >= -fB/2) & (timeVector <= fB/2)) = 1;
psiWavelet = ((pi*fB))^(-0.5).*...                   % Morlet function
exp(2*1i*pi*fC.*timeVector).*exp(-timeVector.^2/fB); % Morlet function
%      FFT plot by matlab bulit-in FFT function
     Nfft =10*2^nextpow2(length(timeVector));
    FFT =fftshift(abs(fft(psiWavelet,Nfft)));
    freqs=[0:Nfft - 1].*(fS/Nfft);
    freqs(freqs >= fS/2) = freqs(freqs >= fS/2) - fS;
    freqs=fftshift(freqs);
    figure(2);
    subplot(1,2,1)
    plot(freqs, FFT);
    xlim([-1  5]);
    xlabel('Frequency / Hz');
    title (sprintf('Fourier Transform'));
%     FFT plot by its direct fourier transfrom function
    f_psi=exp(-(pi^2*fB)*(freqs-fC).^2);
    subplot(1,2,2)
    plot(freqs,f_psi)
    xlim([-1  5]);

If you approximate the Fourier transform

$$X(f)=\mathcal F(x)(f)=\int_{-\infty}^\infty x(t)\,e^{-2\pi j\,ft}\,dt$$

by the discrete Fourier transformation for by sampling on the time segment $[-T,T]$ as

$$X(f_n)\approx \sum_{k=-N}^{N-1} x(k\tau)\,e^{-2\pi j\,f_nk\tau}\,\tau=s[n]\,\tau$$

with $T=N\tau$, $f_n=n/(N\tau)=n/N*f_s=n/T$, $n=-N,...,N-1$, $s$ the result of FFT on the sampled $x$ sequence after shift, ...

then you have to multiply the result of the FFT method by $\tau=T/N=1/f_s$.


Deriving the Resolution Equation of an FMCW RADAR I am interested in the equation that depicts the resolution of an FMCW radar. I write the same below.

$\Delta d = \frac{c}{2b}$

First I would like to have an intuitive explanation to the above derivation.

Secondly, if possible, I would like to know why the chirp's gradient doesn't play a role in determining the resolution of the FMCW radar. Does it not matter how fast I sweep the $BW$? If not so, what impact does the chirp gradient have?

Definitions:

$\tau$: round-trip delay.  $\Delta \tau$ = delay resolution.
$r$: range.  $r=c\tau/2$.  Range resolution $\Delta r = c \Delta \tau / 2$.
$B$: chirp bandwidth
$T$: chirp sweep time
$R$: chirp rate.  $R=B/T$

Basic steps:

We estimate the delay by estimating the frequency shift between the received and transmitted signal.  Typically the resolution of a frequency estimate is the inverse of the observation time, which in this case is $1/T$.  Why?  Think of taking the Fourier Transform of a tone with a rectangular window T.  The window will have nulls at $f-1/T$ and $f+1/T$, where $f$ is the tone's frequency.
Because the chirp rate is $R=B/T$, the time resolution corresponding to a frequency resolution of $\frac{1}{T}$ is $\frac{1/T}{B/T} = \frac{1}{B}$.
Plugging this into the range resolution equation, we have $\Delta r = c/(2B)$.

Intuitively, we improve range resolution by (1) improving resolution of the frequency shift estimate (i.e., increasing sweep time), and (2) decreasing the time delay corresponding to a given frequency shift (i.e., increasing chirp rate).  To do one without sacrificing the other, we need to increase the bandwidth.  If we double the chirp rate but cut the sweep time in half, thus keeping the bandwidth the same, the time delay per frequency shift decreases but the frequency resolution gets worse.  Likewise, if we cut the chirp rate in half but double the sweep time, the time delay per frequency shift increases while the frequency resolution improves.  In either case, the range resolution does not change.


What is the relation between time, frequency, spatial and spectral domains? What is the relation between time, frequency, spatial and spectral domains?

I think that time and spatial domains are the "original" domains (e.g. if the signal is a greyscale image, then the domain would e.g. be all possible grey intensities for all pixels), whereas the frequency and spectral domains are the "Fourier" domains, that is, the domains of the functions after having applied the Fourier transform to the original function or signal. Is this right?

I also think that we use spatial when the input, in the original domain, are images, whereas we use "time" in a more general sense, but maybe we should be using it only when the original domain only represents time (we could say a similar thing regarding the frequency domain). The spectral domain, I think, is a more general term, which refers to the domain after having applied the Fourier transform, so it can e.g. refer to the frequency domain. Is this right?

I am more interested in knowing the meaning of these concepts in the context of image processing and computer vision. I think it would be interesting to have a diagram that explains the transitions between these domains (using e.g. the FT).

A signal (e.g. a time series or an image) is a function $f: A \rightarrow B$, where $A$ and $B$ are respectively the domain and codomain of the signal $f$, which is also often denoted by $x$. A time series signal is often denoted by $x(t)$ to emphasise that the signal $x$ depends on the time $t$.

In general, the expressions "time domain" and "spatial domain" refer to the domain of the original or given signal, which can 1D or 2D (or of higher dimensions). The expression "spatial domain" can be used to more specifically refer to the domain of a given image (which is a subset of the $2D$). On the other hand, "time domain" is a more general term and can be used to refer to the domain of a time series signal or an image, even though an image does not really depend on time (but, by convention, we can still refer to the domain of the image using the expression "time domain").

In general, the expressions "frequency domain" and "spectral domain" refer to the domain of the transformed signal using a Fourier transform. In other words, the domain of the Fourier transform of a given signal, which is another signal, can be referred to as any of these expressions. It can also be called "Fourier domain". So, "frequency domain", "spectral domain" and "Fourier domain" are synonyms. However, the expression "frequency domain" is often used in conjunction with the expression "time domain".


Are Amazon Echo ads embedding an ultrasonic signal to suppress recognition of the Alexa keyword? Related to High audio frequencies in television broadcast

My cell phone samples at 48000 samples/sec so the FFT app I installed can only detect frequencies up to 24 kHz.

If anyone has an audio frequency analyser that can go higher, would you please use it on a digital TV advertisement for Amazon Alexa and see if they are embedding any sort of ultrasonic signal to tell the device not to respond to the keyword?

When it first came out it responded to the word Alexa on TV all the time and now it almost never does.  My guess is an ultrasonic signal is being used to suppress the keyword recognition.  No-one has talked about this online or in the press that I can find, so I'd like some confirmation if this is true.  (I can't see any other way they would be doing it that doesn't involve sending samples to Amazon which they should not be doing since the recognition light is not activating at all for these events)

So my question is whether an ultrasonic signal is present in Alexa TV ads and if so, what is it composed of?

Most speaker/amplifier systems will have analog filters that cut off far below 24 kHz.

In fact, with 44100 kS/s being a very common sampling rate, you can be certain that you can't transport any signal at frequencies above 22.05 kHz; Nyquist doesn't only apply to the digitizing process, but also throughout the rest of the signal processing chain. Since digital spectra are periodic with the sampling rate, any digital-to-analog converter must have so-called reconstruction filters that cut off at half the sampling rate on each channel.

In other case, what your phone can't hear is something that your TV cannot emit, either. (At least not intentionally – analog filters, especially affordable ones, are far from perfect)

My guess that the software simply got better at determining the "intonation" of commands. Luckily, for all Latin- and Germanic-based languages, the imperative is not that hard to discern acoustically – if I had a guess how most of an infant's and a dog's understanding of human language starts, it's that – and you can tell someone commanding "Alexa, do my dishes!" from someone noticing it in the center of a sentence as object like "I can't trust Amazon enough to let Alexa listen in on confidential conversations."

Another relatively common thing is Steganography – in the case of Audio that's usually called Watermarking – the art of hiding info in a signal in a way that is recoverable, but not noticeable to the observer. It's been used for easily one-and-a-half decades now[citation needed] to track bootlegged copies of movie recordings etc.


Simulation of a pressure field developed by a song Greetings!

What I am doing:

In short, I want to write a program that is fed a song as a wav file, and spits out an animated contour plot indicating the relative pressure/loudness/quantity indicating amplitude at various points in a plane. Each speaker is modeled as a point source a set distance from each other hanging in an infinite void. My weapon of choice is Matlab.

How I am doing it:

Poorly. Song is read and split into channels, then cut into little slices. Each slice is FFT'd and the amplitudes and phases are split. For each point in the simulation space every frequency is taken, wavelength calculated and used in combination with the phase to determine the contribution of that frequency to the dimensionless loudness/pressure/sound intensity indicating quantity in that cell (the exact unit isn't that important - explained in next section). Do this for each slice of song, draw contour plots for each slice, and animate!

What my problem is:

The above seems kinda clumsy, and my preliminary results tell me I'm likely neglecting something huge or making some kind of fundamental error. Please note the idea here isn't to make a simulation for speaker setup or any kind of real physical use, I want to draw pretty pictures of what songs do with some basis in reality (specific application: psytrance). This is why the units are not that important, whether I use the normalised voltages in the wav file or convert it to dB or pressure... apart from looking a bit different in lin or log space the relative difference should be the same, at least for my pretty-picture-drawing-purposes. Once the program works as planned some number-fudging will occur to make it look better.

My Questions

1) Can someone suggest a better way of doing this? I feel like I'm doing this in a too complex manner, but I don't even know where to begin looking for a different method to accomplish my goals.

2) Am I right in doing phase corrections with the values given in the FFT? Gut says the error is here - Fourier coefficients are returned as complex numbers a+bi - for the use I have described here, how would I use this to yield "amplitude" and phase angle? It's a simple question... but humour me, self-kicking might follow your answer.

Edit: On second thought this may not be the most appropriate place to post this (at least for part 1 of the question), if someone can suggest a better place, please do.

For a perfect isotropic point source the pressure or air displacement at a point is just the pressure/displacement stimulating signal delayed by the time it took the disturbance to travel to that point.

That means you get spherical shells of constant displacement, and the displacement is the signal at the time $t-d/c$ where $d$ is the distance from the source and $c$ is the propagation speed.

Addendum:

If you do not require the full resolution of the pressure field that would allow you to actually reconstruct the sound from the field then you may want to average the pressure over spatial volumes of constant size. That is not equivalent to a time average of the power, because the influence of spatial coarse graining on the pressure history depends on where you are. Only in the far field where the radius of the constant pressure sphere is much larger than the typical coarse graining scale, time averaging ans spatial averaging become equivalent.

So it depends on what you want exactly and what you need it for. The physics involved is very simple for a perfect point source, but gets very complicated as soon as you have participating objects. Maybe you can give a more detailed description of what you are trying to achieve with this.


Differentiation of discrete sinusoidal signal What I understand from differentiation of discrete sinusoidal signal is $y(n) = \alpha(x(n)-y(n-1))$ where $\alpha$ may be a factor depending on sampling frequency and signal frequency. Now y(n) represents signal which is $90^\circ$ phase shifted sinusoid.

My question is why to go to Hilbert transform and other complex methods to shift a signal $90 ^\circ$. My concern is with $50\,\text{Hz}$ signal only till $40$-th harmonic.

Hilbert transform has the property of harmonic correspondence (Vakman 1996), meaning that the Hilbert transform of a cosine wave is a sine wave etc. So the Hilbert transform is a 90 degree phase shift. In fact the Hilbert transform in the Fourier domain looks quite similar to the differential operator:

$$\mathcal{F}\left[\tfrac{d}{dt} f(t)\right](u) = -iw \mathcal{F}[f(t)](u)$$
$$\mathcal{F}\left[\mathcal{H}[f](t)\right](u) = -i\frac{w}{|w|} \mathcal{F}[f(t)](u)$$

We can see that for a sinusoid of frequency $w_0$, the derivative is the Hilbert transform scaled by $w_0$. e.g. Hilbert transform of $\cos(n t) = -\sin(n t)$, derivative of $\cos(n t) = -n \sin(nt)$. So careful choice of alpha in your formula will give you the phase shift, while compensating for the scaling by frequency, which is the same as the Hilbert transform.

If you have a non-sinuosidal signal, then the 90 degree phase shift must be applied to every sinusoidal component. This is not the same as differentiation, and you couldn't use your formula. e.g. Hilbert transform of square wave (from wikipedia):


What algorithm can I use to "match" or "align" (translate and scale) two vectors of different size in C#? Let's say I have two vectors A and B, of lengths, m, and n respectively. They are not "signals" but a collection of ordered time marks or indexes from another two "original" longer vectors (let's call them X and Y).
I try to compare A and B to "estimate" the similarity between X and Y. If they are aligned I can use a simple Jaccard index to do it. However, X and Y (and so A and B) can be delayed and/or steched with each other. I've assumed that if X and Y are similar, then B should be a linear function of A (B = c*A+d, c and d scalars). Besides, some elements could not match exactly or be missing in the other vector because A and B are of different lengths.
For instance:
A = [1 6 12 25 36 37 40], B = [3 13 25 73 82]
Should match the following pairs:
C = [[1,3],[6,13],[12,25],[36,73],[40,82]] -> B = 2*A+1
What algorithm or method can I use to find the gain c and offset d to get the best alignment between A and B? In other words, the parameters that maximize the amount of common elements between both these vectors.
Typically A,B are around 150 elements, and 0.8 < c < 1.2.

A little context: I'm trying to make a song dedupe application in C#. In fact, X and Y are audio firngerprints (chromaprints) made with AcousticID, which are basically binary matrixes (every row is a 32-bit integer). This questions could be directly solved to align X and Y but I think it could be harder due to the "noisy" nature of the prints (sample) and the larger size of vectors.

The most similar method to your thought is finding a "linear regression". Hence, you can use a c# library or code to compute "linear regression" for the points of in A|B elements pairs.


How can I reconstruct a 3D structure from a set of images? I have a set of contours (set of line segments) grouped in the following way:

$S_i = \{I^0, I^\frac{\pi}{4}, I^\frac{2\pi}{4}, \ldots, I^\frac{7\pi}{4} \}$

where

$S_i$ denote sequence of photos of one concrete object.
$I^j$ denote an image, with $j^{th}$ point of view ($j=0$ means front view).

Here is example of $I^\pi$ (rear view):

How can I reconstruct 3d structure of object with given $S_i$?

Can someone point me to some papers or even give me some keywords? I know that there are a lot of articles that operate with clouds of points and so on, but those don't work as I'm operating with lines.

Actually it is quite a hard topic. Classical multi-view 3d reconstruction deals with point matching in the first place, i.e. find the same point on every image. Given the camera (view) parameters for each image, the original 3d point can be reconstructed. (Using a laser or a projector the scene can be lit so the matching can be done relatively easily.)

The bible of the field is Multiple View Geometry in Computer Vision by Hartley and Zisserman

In the book there is a section about the trifocal tensor, which is a multilinear constraint between 3 views. It contains not just point but line correspondence constraints as well. It can be used for building reconstruction very well.

So your contours should be matched at the first place, and maybe can be reconstructed knowing the camera parameters (camera calibration is also covered in the book). Then you will have contours in 3d but nothing more. For real surfaces you have to do dense point matching. Though the tensor I mentioned look good it is used for straight lines and I am sure that a modern car has curved lines all over.

I do not know how you got those contours but seeing the image you have posted I am quite sceptic about the robustness of that algorithm, so the reconstruction will be poor.

Another method it came to my mind is visual hull or space carving. The contour mathcing should also be done. Running the method on each contour you can have the model.


Measuring group delay through a linear phase filter I'm having a difficult time wrapping my head around a phase/ group delay problem and it corresponding to where you start sampling so I built up a system to test some of these issues. The system passes 2 sinusoids through a lowpass FIR filter with linear phase, then I take an NFFT length fft starting at some arbitrary point of my output data, then measure the phase of each sinusoid. The signal flow is as follows:

[Sinusoids] => [FIR (length 65)] => [FFT] => [phase measurement]

I'm trying to use the fact that group delay can be calculated from the phase information  in the fft as:

but this is where I get lost because it seems that where you start sampling at the output would impact the phase measurement. i.e. you can't start at any arbitrary point. But in a black box system where you can't figure out what the correct output sample to start on is how can you accurately measure group delay? Below is my script. My group delay should be 32, but depending on the variable out_start_idx the group delay changes dramatically. Shouldn't it be that regardless of where you start sampling the output, the phase change between two tones should still be constant?

close all
clear all

fs = 1e6;

%GROUP DELAY of filter (N-1)/2 = 32
%since h = length 65
h = remez(2^6, [0 fs*.05 fs*.15 fs/2]/(fs/2), [1 1 0 0]);

%arbitrarily large fft
NFFT = 2^19;

bin_space  = fs/NFFT;

num_tones = 2;
space = 1000; % spacing in bins
for nn = 1:num_tones
 freqs(nn) = space*(nn)*bin_space;
end

%random starting phases between -pi/2 pi/2
phi_start = rand(1,num_tones)*pi-pi/2;
t = 0:999999;

sig = 0;
for nn = 1:num_tones
 sig = sig + exp(1i*2*pi*freqs(nn)*t/fs + phi_start(nn));
end

out = filter(h,1, [sig zeros(1,length(h)+ 10)]);

%{
figure
subplot(211)
plot((-.5:1/NFFT:.5-1/NFFT)*fs, 20*log10(fftshift(abs(fft(h,NFFT)))))
hold on
plot((-.5:1/NFFT:.5-1/NFFT)*fs, 20*log10(fftshift(abs(fft(sig,NFFT)))))
hold off
subplot(212)
plot((-.5:1/NFFT:.5-1/NFFT)*fs, 20*log10(fftshift(abs(fft(out,NFFT)))))
hold on
plot((-.5:1/NFFT:.5-1/NFFT)*fs, 20*log10(fftshift(abs(fft(h,NFFT)))))
hold off
%}
f = 0:fs/NFFT:fs-fs/NFFT;

% take NFFT amount of samples from arbitrary starting point
% in output
out_start_idx = 1;
data_test = out(out_start_idx:out_start_idx+NFFT-1);

% get index of tones
idx_tones = (space:space:space*(num_tones)) + 1

% take fft of data
fft_out = fft(data_test,NFFT);

% get complex output from each tone bin
tone_IQ = fft_out(idx_tones);

% calculate phase of each tone
phi_tones = angle(tone_IQ)*180/pi;

% calculate group delay
phi_2 = diff(fliplr(phi_tones))
phi_1 = diff(fliplr(phi_start)*180/pi);
tau = (-1/360) * (phi_2 - phi_1)./diff(freqs)
grp_delay = tau*fs

% maybe do some extra processing
phi_out = phi_tones;

% linear fit calculated phases
p = polyfit(freqs,phi_out,1);
y = polyval(p, freqs);

% calculate deviation from linear phase
deviation_from_linear_phase = phi_out -y;

% plot deviation / etc
figure
plot(freqs, y)
hold on
plot(freqs, phi_out, 'r*')

figure
stem(deviation_from_linear_phase)

phi_start = rand(1,num_tones)*pi-pi/2;

The line above illustrates your problem. Phase measurements always need to referenced back to something. An absolute phase measurement without a well defined reference is meaningless.

I believe that the method outlined above that you don't measure the absolute phase but the phase of the output relative to the phase of the input. It probably assumes that you can observe both the input and the output. You can measure at an arbitrary time as long as you measure the input and the output simultaneously .


What is the difference between y[n] and y(k) (using square brackets)? Can I modify y[n] = αy[n-1] + x[n] the same way I would do it with y(k) --> Y(z)?
I need H(z), but I don't know what is the difference between using (k), and (n). I often encountered [k] notation, what is the difference between (k) and [k]?

The notation for discrete-time signals, $x[n]$, was first noticed by me in the original Oppenheim and Schaffer (1975), even though i didn't see that book until 1981.

So before, discrete-time signals in DSP used the same notation in mathematics used for sequences (infinite or not), with the subscript, $x_n$.  In both notations, the argument or subscript play the same role and are only integers.

$$ x[n] \triangleq x_n \qquad \qquad x\in\mathbb{C} \qquad n\in\mathbb{Z} $$

Usually the signal, $x[n]$, is real but it need not be.  But the argument, $n$, must be an integer which makes this no different than a sequence, $x_n$.

With the DTFT and DFT this notation was common pre-O&S,

$$\begin{align}
 X_k &\triangleq \sum\limits_{n=0}^{N-1} x_n e^{-i 2 \pi \frac{nk}{N}} \\
\\
 x_n &= \frac{1}{N} \sum\limits_{k=0}^{N-1} X_k e^{+i 2 \pi \frac{nk}{N}} \\
\end{align}$$

or they would use that awful $W_N^{nk}$ notation.

Now the problem with the "$x_n$" notation is that you would get too many subscripts to deal with, if you had a signal that was a multi-element vector.  And there was nothing in the notation to denote which subscript was the "time-like" variable.  With analog signals, $v(t)$, this wasn't a problem.  If you had a cable with 8 wires (plus ground), you could label them "$v_1(t)$" to "$v_8(t)$" and looking at $v_m(t)$, there was no doubt that it meant the "$m$th element that is a function of time, $t$."  But with "$x_{m,n}$" (as opposed to the O&S way "$x_m[n]$") there is nothing, other than order, differentiating discrete time from the "channel" number.  Add another dimension (such as the 2-dim pixel radiance as a function of time) and you have a worser mess:

$$ x_{r,c}[n] \qquad \text{vs.} \qquad x_{r,c,n} $$

$r$ is row number, $c$ is column number, and $n$ is discrete-time.

So with this notation, there is a parity between a vector of continuous-time (a.k.a. "analog") signals $x_m(t)$ and discrete-time (a.k.a. "digital") signals $x_m[n]$.  So, likely, whatever relationship one might see between the various analog signals $x_m(t)$ in a multi-conductor cable, I might expect to see between the various elements of the digital vector signal $x_m[n]$ and sampling would be no different than with a single-dimensional signal

$$ x_m[n] \triangleq x_m(t)\bigg|_{t=nT} \qquad \qquad T=\frac{1}{f_\mathrm{s}} $$

So this provides for a consistency of notation along with a way to differentiate arguments that must be integers whereas whatever is contained within parenths $f\big( \cdot \big)$ is considered to have no such restriction (could be a real or complex argument).  That is, to me, helpful.

So, if it's a very old DSP paper (pre-1978, perhaps), the notation you would see would almost certainly be $x_n$ for discrete-time signals and the DFT would be $X_k$ for discrete-frequency.  But this has evolved to

$$\begin{align}
 X[k] &\triangleq \sum\limits_{n=0}^{N-1} x[n] e^{-j 2 \pi \frac{nk}{N}} \\
\\
 x[n] &= \frac{1}{N} \sum\limits_{k=0}^{N-1} X[k] e^{+j 2 \pi \frac{nk}{N}} \\
\end{align}$$

and EEs will use "$j$", instead of "$i$" for the imaginary unit.

Math folks might use this notation:

$$\begin{align}
 \hat{f}_k &\triangleq \frac{1}{\sqrt{n}} \sum\limits_{j=0}^{n-1} f_j \, e^{-i 2 \pi \frac{jk}{n}} \\
\\
 f_j &= \frac{1}{\sqrt{n}} \sum\limits_{k=0}^{n-1} \hat{f}_k \, e^{+i 2 \pi \frac{jk}{n}} \\
\end{align}$$


Why does this Fourier series give lower amplitude for max value of signal? I want to approximate below signal using Fourier series on Matlab.

My code is below
close all
clear all
clc
a_0 = 0
t=-10:0.001:10
x_t = 2/5
T = 10;
w0 = 2*pi/T
numOfTerms = 100
plot(t, x_t);

for k=1:numOfTerms-1
   a_k = (1/(w0*i*k*T))*(  exp(-1*i*k*w0)-exp(-3*i*k*w0)+ exp(-7*i*k*w0)- exp(-9*i*k*w0))
   x_t = x_t + a_k * exp(i*k*w0*t)
   plot(t, x_t);
   title(['Number Of Terms = ',num2str(k+1)]);
   pause(0.1);
end

You see that min point of this signal above zero and max point below 1. How can I solve this problem? Why it gives such result?

The formula for the Fourier series is

$$x(t)=\sum_{k=-\infty}^{\infty}c_ke^{jk\omega_0 t}\tag{1}$$

You only use the non-negative indices $k$, but you must also use the negative ones. Since your signal $x(t)$ is real-valued, you can in fact compute everything using only the positive indices, but then you must take two times the real part:

$$x(t)=c_0+2\Re\left\{\sum_{k=1}^{\infty}c_ke^{jk\omega_0 t}\right\}\tag{2}$$

Note that you must add the DC term $c_0$ separately.

This is the corresponding Matlab code:

t=-10:0.001:10;
x_t = 0;
T = 10;
w0 = 2*pi/T;
numOfTerms = 100;

for k=1:numOfTerms-1
   a_k = (1/(w0*i*k*T))*(  exp(-1*i*k*w0)-exp(-3*i*k*w0)+ exp(-7*i*k*w0)- exp(-9*i*k*w0));
   x_t = x_t + a_k * exp(i*k*w0*t);
end

x_t = 2*real(x_t) + 2/5;

plot(t,x_t);


Impulse response of a system in z domain The question 3-23 in the "Discrete-Time Signal Processing - Second Edition" is:

and the solution is:

I cannot understand the solution.  In the second row of the answer when I multiply (-4)  with this denominator:

and then add the result to:

it gives a result different from the numerator of the first row. This multiplication must be equal to

Which is not. Is this solution true? If yes, can anyone help me to understand if?

You are right.

In the solution, the second line which performs the long division

$$ H(z) = -4 + \frac{ 5 + \frac{7}{2} z^{-1} }{1 - \frac{3}{4}z^{-1} + \frac{1}{8} z^{-1} } $$

is wrong and should be corrected as:

$$ H(z) = -4 + \frac{ 5 - 3 z^{-1} }{1 - \frac{3}{4}z^{-1} + \frac{1}{8} z^{-1} } .$$

However, the partial fraction expansion at the following third line is correct; indicating there was a typo in the second line. Finally, the impulse response (which is based on line three) is also correct.


Correlation: A Magnitude or Power quantity? I am confused about the units we would refer to with a correlation result in that if it would be a magnitude or power quantity (and therefore specifically when using ratios in dB would we use 10Log or 20Log?).
I have since confirmed as detailed in this post that the normalized correlation coefficient $\rho$ is a magnitude quantity in that the relationship between SNR and $\rho$ repeated here is:
$$\text{SNR} = 10\log_{10}\frac{\rho^2}{1-\rho^2}$$
Where in this case for an SNR relationship to $\rho$, the correlation involved is the correlation of $x(t)$ being a reference waveform as signal to $y(t)$ which is the same signal with added noise. (So that the SNR of $y(t)$ is determined in this case).
However I also understand and see directly from the math that the autocorrelation (at $\tau=0$ if we are referring to the autocorrelation function) is the variance (scaled by the number of samples), which is a power quantity, and we see the sum of products in the general expression for correlation suggesting a sum of powers.  How do we reconcile all this?
Is it that the normalized coefficient through it's normalization process:
$$\rho = \frac {\operatorname*{cov}(x,y)}{\operatorname*{stddev}(x)\operatorname*{stddev}(y)} $$
Has this converted the power quantity of the numerator back to a magnitude quantity, and therefore correlation on its own (the numerator) IS a power quantity, while the correlation coefficient IS a magnitude quantity? If so, I don't yet quite see how. Or is it something else entirely?

There are different ways to look at this.
In general, it's useful to think of the correlation as a power quantity.

If you slog through the units, the correlation will have power-like units: $V^2$ for example
The Fourier Transform of the autocorrelation is the Power Spectral Density, which is clearly a power-like quantity

So what's happening in your example? Let's assume that we have $y = x + n$ where $x$ and $n$ are mean-free, uncorrelated and have the variances $P_x = \sigma_x^2 = <x^2>$ and $P_n = \sigma_n^2 = <n^2>$. Since they are uncorrelated we have $P_y = P_x+P_n $
The SNR would be $$SNR = 10log_{10}\frac{P_x }{P_n}$$
The argument to the log is a ratio of powers.
The covariance is
$cov(x,y) = <x\cdot (x+n)> = P_x$ which has power-like units. The normalized correlation coefficient comes to be.
$$\rho = \frac{cov(x,y)}{\sqrt{{P_x \cdot P_y }} = \frac{P_x}{\sqrt{P_x \cdot (P_x + P_n) }} $$
$$\rho = \frac{cov(x,y)}{\sqrt{P_x \cdot P_y }} = \frac{P_x}{\sqrt{P_x \cdot (P_x + P_n) }} $$
So that's also a ratio of powers: cross-power to geometric mean of the signal powers. Since it's a ratio it's unitless. Let's square it:
$$\rho^2 =\frac{P_x^2}{P_x \cdot (P_x + P_n) } = \frac{P_x}{P_x + P_n } $$
The most reasonable physical interpretation of this is STILL a ratio powers: signal power to signal power plus noise power. Since the ratio is unitless, the square of ratio is still unitless, so there is no direct conflict here.


Audio Activity Detection and Classification I am starting a new project. Actually, my real intention is to learn speech recognition but for warming up I want to improve in audio signal processing.
In my project, I aim to record sounds except for human speech in my house when I go out. I am planning to divide my project two sections: Audio Activity Detection and Audio Classification. However, some questions appear in my mind.

What is the best techniques for audio classification? Could you supply some informative articles for me? I saw that MFCC and CNN are generally used for speech recognition. They are also valid for audio classification? Or, should I seek another methods? My classification includes baby cry, door knocks, door opening, door bell, ringing home phone, smoke detector etc. Even though their frequencies, in contrast to frequency of human speech, cannot be limited in a specific range, should I study on frequency for classification?

When I say Audio Detection Activity, I imply that sound recording starts if silence is broken, i.e. there is a sound. For preprocessing, I will filter out the noises. Therefore, noise will be part of the silence.:)  I searched few articles and I encountered the topic of Voice Activity Detection (VAD) but it is used to understand whether human speech is present or not in the environment. I am not interested in voices. Instead, I want to collect sounds mentioned before. Is Audio Detection Activity a valid term?

I am waiting for your advices and sources.

Thank you.

Here is an recent introduction video on the subject (2019): Audio Classification with Machine Learning


How can both Gaussian and salt & pepper noise be removed from an image? An image is corrupted by Gaussian noise,then salt & pepper noise. How can you eliminate it and why?

The easiest approach: First apply Median filter for salt and pepper noise, and then use Gaussian blur to eliminate Gaussian noise.


Find output of digital filter given input and transfer function? Hey guys so I have an input sequence that are real values that represent intensity of a column of an image. I also have the transfer function in the Z domain of a 2nd order high pass filter. I have the following transfer function using the bi-linear transform estimate.

Say my input signal is defined as x[n] = 27, 588, 1013, 255, 38. How can I determine the output y[n]?

There are two options, one being I determine X[z] then multiply it by H[z] and take the inverse transform.

The other being I transform H[z] into h[n] and then I convolve it with x[n] using the following scheme .

How would I implement either of these solutions in code (eg. coding an automatic way to generate X[z] from an arbitrary x[n] OR programming a convolution function that implements FOIL for the second option)?

Is one easier than the other? Keep in mind that my signal is 512 samples in length as well, so convolution may be very expensive.

Standard difference equation of a biquad: y[n] = b0*x[n]+b1*x[n-1]+b2*x[n-2]-a1*y[n-1]-a2*y[n-2];

Initialize all x[n], y[n]  to 0 for n < -1

b0,b1,b2 are 0.86 -1.72 0.86 (numerator)

a0,a1,a2 are 1 -1.71 0.73 (denominator)


SNR contains $\log N$ term where $N$ stands for number of samples Sampling the signal$ $N times increases the signal energy by a factor of $N^2$ and the noise energy by a factor of $N$. Why? This explanation is written for
              SNR(dB) = signal peak(dB) – noise floor(dB)- $10\log N$

In the additive model $y=s+n$, when the signal is deterministic, it adds coherently over the "realizations". Hence, its variance $V(\sum s_n) = V(N s) = N^2  V( s)$. And when the noise $w$ is independent identically distributed (IID), then $V(\sum w_n) = NV( n) $. This is a classical result on the  Variance of Uncorrelated Variables.

You can find this classical result detailed in Deriving the SNR for averaged signals. Warning, this does reach a limit when data is quantized, in other words the noise floor does not follow $1\sqrt{N}$, but plateaues.


Is possible reach the DFT if I have the DTFT? My teacher told me that DFT is DTFT sampled, i.e.:

$$X[k] = X(e^{j \omega})\Bigg|_{\omega = \frac{2\pi k}{N}}$$

But, if I have the sine

$$ x[n] = \sin(\omega_0 n) $$

the DTFT is:

$$X(e^{j \omega}) = \frac{\pi}{j} \big(\delta(\omega - \omega_0) - \delta(\omega + \omega_0)\big)$$

(periodic in $\omega$ with period $2\pi$). How is it possible to reconcile this with the DFT for $N$ values for this sine or any signal?

Thanks people. Sorry for the English.

The DFT is a sampled version of the DTFT only for finite length signals. Otherwise, there is no point in comparing the DTFT with the DFT because you can only compute the DFT for finite length (or periodic) signals. Your example $x[n]=\sin(\omega_0n)$ is an infinitely long sequence. For specific choices of $\omega_0$ it is periodic, but nevertheless, you can't obtain its DFT by sampling its DTFT. The operation of sampling a Dirac impulse is undefined because $\delta(\omega)$ cannot be evaluated at $\omega=0$ since the Dirac impulse is not an ordinary function but a distribution.


Automatic detection of noise-only segments in audio This is my first question for this community. My connection to DSP comes via machine learning/deep learning, and I am working in the Python ecosystem.
As a preprocessing step to audio classification, I am looking for a simple method to reduce noise.
Promising noise reduction methods such as 1 expect a noise-only sample as input, then are able to subtract noise. (The approach is also described in an answer of 2)
Getting a noise sample manually is doable for denoising a few recordings. However, it does not scale to processing thousands of recordings with different noise types and volumes.
Can you recommend a method for automatically obtaining a noise sample from every recording? Alternatively, denoising methods that do not require a noise sample as input?

The way I read it, the algorithm works by figuring out where the signal part lives in the frequency domain by using the statistics of the noise part. This means you don't need a noise-only signal for each signal you want to process as long as the statistics of the noise are the same for each signal realization.
For example, if you had 1000 signal realizations that you want to perform this noise reduction on, you would use a noise-only signal to compute the statistics of the noise. Then using those statistics, the algorithm will figure out signal versus part by comparing to a threshold (average noise power or something I would guess).


Are complex exponentials the only eigenfunctions of LTI systems? Is there an example of an eigenfunction of a linear time invariant (LTI) system that is not a complex exponential? Justin Romberg's Eigenfunctions of LTI Systems says such eigenfuctions do exist, but I am not able to find one.

All eigenfunctions of an LTI system can be described in terms of complex exponentials, and complex exponentials form a complete basis of the signal space. However, if you have a system that is degenerate, meaning you have eigensubspaces of dimension >1, then the eigenvectors to the corresponding eigenvalue are all linear combination of vectors from the subspace. And linear combinations of complex exponentials of different frequencies are not complex exponentials anymore.

Very simple example: The identity operator 1 as an LTI system has the whole signal space as eigensubspace with eigenvalue 1. That implies ALL functions are eigenfunctions.


From a physics perspective, why does D/A quantization error result in a noise floor? For the last week or so I have been trying to understand how quantization error results in the noise floor outside of a mathematical perspective and I haven't really had any luck finding a source that discussed quantization noise without using equations to show where quantization error comes from.  I certainly have nothing against using math proofs but I was hoping to understand why quantization noise existed rather than how to calculate it.  I think I’ve been able to work through it but I was hoping someone would be able to double check my work or point me in the right direction if I’m missing something.

So here is how I figure it so far:
Because the sample rate is allowing us to have a perfect recreation of the frequency of the wave the bit depth only determines the range of values available for the amplitude of a wave being converted to a digital signal.

When the original signal gets rounded to the nearest bit it creates slight changes in the amplitude which result in a slightly altered waveform at each sample.  This new waveform is basically the same as the sum of our original signal and a new frequency.  So if we were to duplicate this in an analog signal we could do it by adding additional frequencies at very small amplitudes to our signal 44.1 thousand times a second (or whatever the sample rate is).  The lower the bit depth the bigger these changes to our original waveform will be, and so the larger the amplitude of the added waves will be raising the “noise” floor.

If you sampled a waveform that allowed for a slower sample rate, say sampling a 10 Hz wave 20 times a second.  Instead white noise as a result of quantization error would you hear a quick succession of random frequencies paired with your original wave (which would obviously be to low to hear)?

If you can follow my logic up to this point I have one additional question because in my tests changing the phase of my sampled analog wave doesn't seem to change the frequencies of the quantization error like I thought it would. Am I making a mistake in my test?

For the last week or so I have been trying to understand how quantization error results in the noise floor outside of a mathematical perspective and I haven't really had any luck finding a source that discussed quantization noise without using equations to show where quantization error comes from.

That's because it's a purely mathematical effect, there's no physics involved.

It's technically quantization distortion, not noise, but we simplify the analysis by treating it as a noise source.

But it's really not; it's nonlinear distortion, so the spectrum of the quantization error varies with the signal.  Maybe that's where you're confused?

For instance, if you sample a sine wave that's a sub-multiple of the sampling rate and quantize it:

you will get harmonic distortion only at specific frequencies:

If it's not an exact sub-multiple, the infinite number of harmonics produced by the distortion will alias and produce a more noise-like residual:

Since this sort of behavior is more common with arbitrary signals, it's usually treated like a noise source.  In addition, we usually randomize this distortion by adding a small amount of dither noise to the analog signal before quantization, which flattens out the spectrum, reducing the peak value of the distortion and making the noise unrelated to the signal, so we are actually forcing it to be a true noise source rather than a distortion:

from __future__ import division
from numpy import linspace, round, abs, sin, log10, pi
from numpy.fft import rfft
from numpy.random import randn
from scipy.signal import hamming as window
import matplotlib.pyplot as plt

fs = 5000  # Hz
t = linspace(0, 1, fs, endpoint=False)
f = 25  # Hz
x = sin(2*pi*f*t)
plt.figure(1)
plt.plot(x)
plt.figure(2)

x += randn(fs)/20  # Dither
x = round(x*5)/5  # Quantization

plt.plot(20*log10(1/fs * abs(rfft(window(fs)*x))))
plt.ylim(-100, 0)
plt.figure(1)
plt.plot(x)
plt.margins(0, 0.1)


How to get a Rician fading channel matrix from a Rayleigh fading channel matrix? I have a Rayleigh fading channel matrix H. Let's say the dimension is 2000*128. Is there any way to convert this matrix into a Rician fading channel matrix for example κ = 5, 10, 15, etc.? Is there a formula for that?

A Rician fading  channel is just a Rayleigh fading channel with a specular (that is, unfaded) signal component.  Depending on what exactly you mean by the  Rayleigh fading channel matrix $H$, the answer could be as simple as changing the distribution of each element of $H$ from a complex zero-mean proper Gaussian random variable to a complex Gaussian random variable with mean of the form $A_{m,n} \exp(j\theta_{m,n})$.


Zero-Padding seems to CHANGE FFT of signal! I am having trouble with zero-padding. As I understand it, in theory, zero-padding in the time-domain should 'sinc-interpolate' the FFT of the signal in the frequency domain.

This would mean, that 1) none of the values of the original signal change. 2) the overall profile of the FFT doesnt change- aka the peaks remain around the same point.

However, in simulations this does not work.

For example, I FFT the following sequence

  x=[0 1 0 1 0 1 0 1]

at an imaginary same frequency, Fs=1Hz. With no zero-padding (FFT-8),

  |X(f=.3125)|=4.

With 2x zero-padding (FFT-14),

  |X(f=.3125)|=2.6.

As $N \rightarrow \infty$, $|X(f=.3125)| \rightarrow 0\ $. Also, the peak of the FFT shifts from .3125 to .25 as $N \rightarrow \infty$.

Here is my code. The relavent parameter is ZP, which controls the amount of zero-padding (N=8*ZP):

ZP=0;       % of times of original length of signal to ZP with

x=[0 1 0 1 0 1 0 1];
x=[x,zeros(1,length(x)*2^ZP-length(x))];
X=fft(x);

Fs=1;
L=length(x);
df=.5*Fs/L;
t=(1:L)*(1/Fs);
f=(1:L)*df;

figure(12),
subplot(211), stem(x);
subplot(212), stem(f,abs(X)); xlim([-.02,.52])

test_pnt=.3125;
I=find(f==test_pnt);
title(['abs[X(' mat2str(test_pnt) ')]=' mat2str(abs(X(I)))]);

Your problem is these two lines:

t=(1:L)*(1/Fs);
f=(1:L)*df;

They should read:

t=(0:L-1)*(1/Fs);
f=(0:L-1)*df;

If I modify your script as below (for scilab, not matlab), I get the following plot.  The peaks all line up correctly.

And the peak of the FFT is at 0.25 not 0.3125.

for ZP=0:5       // of times of original length of signal to ZP with

    x=[0 1 0 1 0 1 0 1];
    x=[x,zeros(1,length(x)*2^ZP-length(x))];
    X=fft(x);

    Fs=1;
    L=length(x);
    df=.5*Fs/L;
    t=(0:L-1)*(1/Fs);
    f=(0:L-1)*df;

    figure;
    subplot(211), plot(x);
    subplot(212), plot(f,abs(X)); //xlim([-.02,.52])

    test_pnt=.3125;
    I=find(f==test_pnt);
    title('abs[X(' + string(test_pnt) +']='+ string(abs(X(I))));

end;


What should be the output of the pdf for multivariate normal distrubition? Let's think we have a multivariate normal distribution with
$$\mathbf x_1=\begin{bmatrix}x_{11}\\x_{12}\end{bmatrix},\quad \mathbf x_2=\begin{bmatrix}x_{21}\\x_{22}\end{bmatrix},\quad\mathbf m=\begin{bmatrix}m_1\\m_2\end{bmatrix}, \quad\text{and}\quad\mathbf S.$$
Where $\mathbf m$ is the mean, and $\mathbf S$ is a $2\times 2$ covariance matrix.

What should the result of the PDF be? Should it be a $2\times 1$ array or a single value?
Can  any of its values be bigger than $1$?

The PDF $p$ of a, say, bivariate Normal distribution with mean $\mu$ and covariance $S$ is a function that takes two parameters and has a scalar-valued output:

$$ p: \mathbb{R}^2 \mapsto \mathbb{R}$$

and is given by

$$p(x_1,x_2) = \frac{1}{2\pi |S|}\exp\left(-\frac{1}{2}(\vec{x}-\vec{\mu})^TS^{-1}(\vec{x}-\vec{\mu})\right)$$

where $\vec{\mu}\in\mathbb{R}^2$ is the mean of the distribution and $S\in\mathbb{R}^{2\times 2}$ is its covariance matrix. $\vec{x}=[x_1, x_2]^T$ is the function input and $p(x_1, x_2)$ is the function output.

So, the output is a scalar (a single value). Furthermore, the PDF can become larger than one for some input $\vec{x}$.


Analytical Expression for Convolution of Two 2D DFTs I am trying to do an image analysis problem on some images, and I want to calculate some things in the frequency domain which are giving me problems. Say I have an (discrete) image $I(x,y)$ of size $N_x\times N_y$ (in practice this size is about $150\times40$ pixels). I want to calculate

$$DFT\left[\frac{\partial I}{\partial x}\cdot\frac{\partial I}{\partial y}\cdot\frac{\partial^2 I}{\partial x \partial y}\right]$$

and express this in terms of the DFT of $I$, which is an $N_x\times N_y$ matrix with entries $\hat{I}_{k,\ell}$. The product is to be interpreted pointwise. I discretize the derivatives with central differences, e.g.

$$\frac{\partial I}{\partial x} = \frac{I(x+1,y) - I(x-1,y)}{2\Delta x}$$
and so I can compute the DFT of various derivatives using the shift properties of DFT. I'm doing this in the frequency domain because other parts of my problem involve 4th derivatives which I think will be way too sensitive to compute in real-space, and I hope the frequency-space calculations will be more stable/robust.

I've attempted to calculate a slightly simpler case, involving only one derivative. My attempted solution is to use the convolution duality properties of the DFT and express multiplication in real space as convolution in frequency space. My attempted expression is below:

$$DFT\left[\left(\frac{\partial I}{\partial x}\right)^2\right]_{k,\ell} = DFT\left[\frac{\partial I}{\partial x}\right]*DFT\left[\frac{\partial I}{\partial x}\right] = \frac{-1}{N_x(\Delta x)^2}\times\sum_{m=0}^{N_x-1}\sin\left(\frac{2\pi m}{N_y}\right)\sin\left(\frac{2\pi (k-m)}{N_y}\right)\hat{I}_{m,\ell}\hat{I}_{(k-m),\ell}$$

However, I'm not sure if I'm properly generalizing circular convolution to more than one dimension, and even if I am, I don't know how to interpret the convolution property with a mix of derivatives (e.g. $DFT[\frac{\partial I}{\partial x}\frac{\partial I}{\partial y}]$). What am I summing to, $N_x-1$ or $N_y-1$? I suspect I might have something wrong in generalizing circular convolution from 1 to 2 dimensions.

  In general, how should I interpret circular convolution in 2D? And how
  can that shed light on calculating this particular DFT?

Convert your PSF to OTF. Then, multiply your OTF by the image (in which it should be in the frequency domain as well as your OTF).


Dimension error using findpeaks function in Matlab on Power Spectral Density to calculate SNR of audio data I have a set of audio files for which I would like to determine a SNR metric. I tried to implement in Matlab as follows:
% Calculate the FFT and PSD of each of the signals in the signal vector.
afftvectors = {};
hfftvectors = {};
psdvectors = {};
halffreqvectors = {};
powdbvectors = {};
signalpeaks = {};
signalpeaklocs = {};
noisepsd = {};

for k = 1:numel(signalvectors)
    afy = abs(fft(signalvectors{k}));
    afftvectors{k} = afy;
    % Due to the symmetric nature of FFT, only half the spectrum is actually required.
    hfy = afy(1:length(afy)/2+1);
    hfftvectors{k} = hfy;

    py = (1/(freqarray{k} * length(signalvectors{k}))) * abs(hfftvectors{k}).^2;
    % since we take only half of the spectrum, we can conserve the energy by multiplying
    % by 2. Zero frequency (DC) (index 1) and the Nyquist frequency (last index, i.e, end)
    % do not occur twice so, we'll skip the multiplication operation on them.
    py(2:end-1) = 2 * py(2:end-1);
    psdvectors{k} = py;
    halffreqvectors{k} = 0:freqarray{k}/length(signalvectors{k}):freqarray{k}/2;
    powdbvectors{k} = 10 * log10(psdvectors{k});

    % Find peaks
    [signalpeaks{k},signalpeaklocs{k}] = findpeaks(psdvectors{k});
    noisepsd{k} = psdvectors{k} - signalpeaks{k};
end

So, the idea was to get the peaks of Power Spectral Density vector. I hit a dead end when I got the output where 'signalpeaks{1}' is a vector of size 13235x1 but my 'psdvectors{1}' itself is a vector of size 36001x1. All the elements in cell array 'signalvectors' has a sample size of 72000. What exactly is going wrong here? Am I using the wrong function in findpeaks() or something else needs to be done?

Firstly, here is how you can fix what you are tryin to do:
noisepsd{k} = psdvectors{k};
noisepsd{k}(signalpeaklocs{k}) = noisepsd{k}(signalpeaklocs{k}) - signalpeaks{k};

Secondly, here is why I think it is not the greatest idea: You are finding the peaks and subtracting them from the spectrum. This is actually equivalent to setting the points where you have local peaks to zero. In effect, you will have a spectrum that is equal to zero at all points where you have local peaks and equal to signal plus noise at all other points. I fail to see why this could be a good estimate of the noise power spectral density.
Thirdly, if you ask me how to do better, well... finding the SNR of a signal is like asking how to get A and B from C = A+B. You cannot! The problem is not unique. If we have a noisy signal, how would we know what is signal and what is noise, how to even properly define "noise" here? In general you cannot, because both could be stochastic. You can only find an approximation and this approximation will only work well if you know something about your signal. For instance, if you had a pure tone in noise, the problem becomes solvable. We can estimate the tone parameters (frequency, amplitude, phase) form the signal, remove it and treat the remaining part as noise. It works in this case because we know how to differentiate signal and noise, the signal should be tone-like. However, in general audio signals need not be tone-like at all.
The only answer I can give is studying your signals closely and trying to come up with some heuristic that describes them well. Maybe your signals are bandlimited, i.e., you have some regions in the frequency range where you expect no signals. Then you could use those to estimate the noise floor. You would need to assume thought that the noise floor is constant (i.e., the noise is white).
All this will only work if you know a thing or two about your signals. If you know nothing it's back to guessing A and B from C = A+B.


FFT function of Gaussian pulse I need gaussian pulse fft. I set the sigma of gaussian pulse 1e-10 and define the gaussian function in time and use the fft. as you know the amplitude of should be unit in frequency domain. but my result shows that the amplitude of the pulse in frequency domain equal to 1/end_time. for example, if i define the time as t=0:1/fs:1e-07 then the amplitude of X(f) is 1e7 and I changed the end_time then the amplitude respectively change.

this is my code :

clear all;
clc;
fs=1e10; %sampling frequency
sigma=1e-10;

end_time=1e-7;

t=0:1/fs:1e-7; %time base

t_s=1e-8;
variance=sigma^2;
x=1/(sqrt(2*pi*variance))*(exp(-(t-t_s).^2/(2*variance)));
subplot(2,1,1)
plot(t,x,'b');
title(['Gaussian Pulse \sigma=', num2str(sigma),'s']);
xlabel('Time(s)');
ylabel('Amplitude');

L=length(x);

X =fftshift(fft(x,L));
Pxx=X.*conj(X)/(L*L); %computing power with proper scaling
f = fs*(0:L-1)/L; %Frequency Vector

subplot(2,1,2)
plot(f,abs(X)/(L),'r');
title('Magnitude of FFT');
xlabel('Frequency (Hz)')
ylabel('Magnitude |X(f)|');
xlim([0 1e10])

The FFT algorithm is not usually normalized while giving the output. Hence it is not unity.It gives as
The Fourier transform is taken as
$$
 F(e^{jw})=\sum x[n] e^{-jwn}
$$
Thus you could see the sum of x's when w=0 case instead of normalized one


Are there any open source biomedical signal processing projects? I'm in a biomedical engineering master's program and am mostly interested in signal processing. Unfortunately, I feel that I'd like to get deeper into the subject, and my university can't accommodate that request -- there's a limited number of courses, and the professor I've asked doesn't know of any local project that I could join.

Are there any open source projects related to biomedical/bioelectric signal processing?

If not, could you recommend any other type of beneficial extracurricular activities? I could probably spend 1-2h a day on this project.

For image processing in biology, there is the FIJI project, a very worthwhile endeavor. They have a list of project ideas as well, maybe there's something that's perfect for you?


anti-spoofing systems I am working with face anti-spoofing. The database consists of a number of real face images and a number of fake images. After pre-processing and feature extraction, I have a feature vector for each image. For classification, the Euclidean distance measure is used. The test image vector is compared with all real image vectors stored in the database, and the minimum distance is chosen (min.1) . The same vector is compared with all fake image vectors and the minimum value is chosen  (min. 2). The minimum rule is applied again to choose the minimum of both min.1 and min.2 to identify the correct class of the test image. My problem is in plotting the ROC curve.  To draw a ROC curve, only the true positive rate (TPR) and false positive rate (FPR) are needed (as functions of some classifier parameter). The TPR defines how many correct positive results occur among all positive samples available during the test. FPR, on the other hand, defines how many incorrect positive results occur among all negative samples available during the test.
I do not have a threshold which can be changed, how can I plot this curve?

I am not really sure what is meant by "fake face" but for the purposes of this response, let's consider that you have two classes "REAL" and "FAKE" and representative example data points from each class are positioned in some $n$ dimensional space, with $n$ denoting the number of features.

In the described situation, the Minimum Distance Classifier (MDC) is used to discover the decision boundary between the two classes (the MDC can work with more than two classes if required but let's just constrain it to two here).

The MDC discovers the decision boundary by first estimating the centroid of each class by averaging the feature values. Once this is done, the classifier is essentially "trained". The knowledge of what it means to belong to one or the other class is basically depending on the position of the centroids in multidimensional space.

The next step is to find the line that passes through the two centroids and its vertical that is further at the half-way point between the two centroids. This vertical line at the half-way point is the actual decision boundary. Anything landing to the right of this line is classified as belonging to one class and anything landing on the left of this line (or on the line itself) is classified as belonging to the the other class.

For the purposes of illustration, I will be re-using one of the figures from a previous link here, let's call the "O"s -> REAL and the "X"s -> FAKE:

So, here is a simple "exercise". Forget about the vertical line being at the HALF-WAY point between the two classes and let's try and draw it a bit to the LEFT, that is, closer to the REAL class.

Obviously, this will classify some of the examples of Class REAL as Class REAL (True Positives), some of the examples of Class REAL as Class FAKE (False Negatives), some of the examples of Class FAKE as Class REAL (False Positives) and some of the examples of class FAKE as Class FAKE (True Negatives).

This "exercise" just gave us the 4 numbers we need to calculate the two metrics of False Positive Rate (FPR), True Positive Rate (TPR) that we need to derive one point of the ROC curve at some threshold value.

So, in the case of the MD classifier, your threshold is equivalent to the position of the decision boundary between the two classes. By "sweeping" that decision boundary all the way from the centroid of one class to the centroid of the other class and deriving FPR, TPR for each value, you will be able to derive the complete ROC curve.

Another way to "visualise" this process with respect to the textbook examples of 'sweeping the threshold value', is to imagine that you rotate your feature space clockwise, until the decision boundary is vertical to the x-axis. If you were to do this, you would be practically eliminating one of the features and you would be "back" to the one dimensional case that is depicted here (for example), with each "gaussian" curve representing the density of the points above that x-axis value. (This is suggested here just as an aid, in case it helps you to make the connection between the position of the decision boundary and the threshold. You don't have to actually rotate the feature space).

Hope this helps.


best similarity measure for audio classification I'm doing a project on audio classification, mainly focused on Bird sounds.

I would like to know which is the best similarity measure that I can adopt to check test signal with training data.

I intent to do the project in matlab.

I would like to know which is the best similarity measure that I can adopt to check test signal with training data.

Welcome to engineering. There is no universally "best" approach, as it always depends on your specific models.

Mathematically, as soon as you can define what "similarity" should mean, you might be able to derive something.

Assuming that bird's auditory senses kind of work similar to human hearing, I'd try with the usual audio comparison models, probably a basic MFCC thing.

As a comment: A former advisor of me once said:

  A week in the lab can easily save you up to two hours of reading literature.

Really, there's a whole lot of papers, textbook introductions and wikipedia articles on sound comparison. Start with anything, and read forward. That's going to be faster than asking for a single answer to your acute problem without getting the bigger picture of how things work.


How to construct Polyphase (Frank) codes in MATLAB I have written a piece of code in MATLAB to generate Frank codes - Pastebin link for pulse compression.

I know that Frank codes are stepped frequency approximations to Linear Frequency modulation(LFM).

LFM Parameters are obtained here:

PW = 100e-6;% Pulse Width
PCR = 100; % Pulse Compression Ratio
CPW = PW/PCR;% Compressed pulse width
Fm = 2/CPW; %Frequency sweep

Stepped frequency approximation to LFM:steppedFreq = 0:Fm/(N-1):Fm;
The phase matrix is generated in the line:Matrix = mod(Matrix*incPhi,360);
Frank code is generated inside the for loop: FCP = [FCP,exp(1j*(2*pift + Matrix(i,p)*pi/180))];
The autocorrelation in dB is calculated here:AutoCorrelationdB = 20*log10(abs(AutoCorrelation)/max(abs(AutoCorrelation)));

My question:

The Peak side lobe level obtained is only around $3.5\textrm{ dB}$. Am I following a correct method for generating Frank codes? The peak side lobe level is given by the formula 20*log10(1/(N*sin(pi/M))) as per the book Principles of Modern Radar: Basic Principles by M.A. Richards - Page 825 which is equal to $29\textrm{ dB}$ for my case($N = 100,M = 10$).

% created by MUSTAFA SAMI AHMED
% Email mustafa_sami87@yahoo.com
clc;
clear all;
M =100;
L=sqrt(M);
bv = 0:L-1;
Matrix = bv'*bv ;
anglout = radtodeg(pi);
% alfa= (4*pi)/(L);

% Matrix_new = alfa*Matrix;
incPhi = (2*pi)/L; %incremental phase change
Matrix = mod(Matrix*incPhi,2*pi);

b=reshape(Matrix,1,[]); % matrix 1D size 1xM
for p = 1:M
  d(p) = exp(1j*b(p)) ;
end

AutoCorrelation = xcorr(d);

AutoCorrelationdB =10*log10(abs(AutoCorrelation)/max(abs(AutoCorrelation)));
 figure(1);
 stem(b);
%   x=wden(AutoCorrelationdB,'heursure','h','mln',10,'sym10');
  %x=wiener2(AutoCorrelationdB, [3 80 ]);
  figure(2);
plot((AutoCorrelationdB));


An Interesting Model with Unknown Orthogonal Design Matrix Suppose the   multivariate one-way anova  model for the raw data , i.e.
$$
\label{Example_model_1}
\mathbf{y}_{ij}=\mathbf{\mu}+\mathbf{z}_i+\mathbf{e}_{ij}, ~~ i=1,\ldots,m,~~j=1,\ldots,n_i,~~~~~~~~~~~~~~~~(1)
$$
where  $\mathbf{\mu}$ is the unknown overall mean of $\mathbf{y}_{ij}$,
$\mathbf{z}_i$ is the unknown random effect at the $ith$ level,  $\mathbf{e}_{ij}$ are the errors,  and
$$
\mathbf{e}_{ij}\mid \mathbf{\Lambda}_0 \overset{iid}{\sim} N_k(\mathbf{0},\mathbf{\Lambda}_0), ~~
\mathbf{z}_i \mid \mathbf{\Lambda}_1 \overset{iid}{\sim}  N_k(\mathbf{0},\mathbf{\Lambda}_1).
$$
Howerever, the raw data is not observed; instead $\mathbf{y}_{ij}^* = \mathbf{\Gamma} \mathbf{y}$ is observed, where $\mathbf{\Gamma}$
is an unknown orthogonal matrix (i.e., the vector observations have been randomly rotated, with
a common unknown rotation).
So the model (1) becomes
$$
\label{model_1_1}
\mathbf{y}^{*}_{ij}=\mathbf{\mu}^*+\mathbf{z}_i^{*}+\mathbf{e}_{ij}^{*}, ~~ i=1,\ldots,m,~~j=1,\ldots,n_i,~~~~~~~~~~~~~~(2)
$$
where $\mathbf{\mu}^*=\mathbf{\Gamma}\mathbf{\mu},$ $\mathbf{e}_{ij}^*\mid (\mathbf{\Lambda}_0,\mathbf{\Gamma})
\sim N_k(\mathbf{0},\mathbf{\Gamma}\mathbf{\Lambda}_0\mathbf{\Gamma}'),$ and
$\mathbf{z}_i^* \mid (\mathbf{\Lambda}_1,\mathbf{\Gamma}) \sim N_k(\mathbf{0},\mathbf{\Gamma}\mathbf{\Lambda}_1\mathbf{\Gamma}'). $

Our aim is to estimate $(\mathbf{\mu},\mathbf{\Gamma},\mathbf{\Lambda}_0,\mathbf{\Lambda}_1).$

Question:
I can't find the relative background of the above model. Please give me some reference or suggest.

We find such signal models in radar literature. For example, $\Gamma$ represents the transmitted signal information matrix, $\mu$ is the deterministic unknown vector representing the reflections from the target, $z_i$ is a random vector generated by the clutter surrounding the target, and $e_i$ is the receiver noise. In statistical literature, such models are called linear mixed models. You can make several assumptions, such as orthogonal signal matrix, unknown clutter covariance matrix etc., and solve a hypothesis testing problem where the goal of the detector is to detect if there is a target or not. In other words, $\mu = 0$ vs $\mu \ne 0$.


Why does it seem most of people will optimize the downlink rate,not the uplink rate? I search for some papers about energy harvest,SWIPT and optimization.And i found that there is just one paper to optimize the uplink rate,the others are for optimizing the harvested power,power budget in the transmitter or downlink rate.

So i want to ask why ? why do most of researchers won't optimize the uplink rate?is there any reason or difficulty?

Because in IoT,we may control lots of things,and of course we do need to know some information of the devices or what happened with the devices,that is ,we do need to know some feedback about the device,so we need  the device uplink the information back to the transmitter,so i think optimizing the uplink rate is still a worth research topic,why  why do most of reseachers optimize the downlink rate instead of uplink rate?

By the way,i also want to ask is the throughput the same as the rate?what is the difference between the throughput and capacity?

This simply reflects the load asymmetry of real world applications. On most networks, traffic is dominated by content streaming, which is a download activity. Netflix alone accounts for 15% of the world's entire Internet traffic.


Why optical flow ? why not tracking I have a question that might be stupide !

Supposing that I try to detect object moving in a video or human action.

Many works are based on optical flow computation.

My question is why using OF is the best for motion detection. A tracking of pixels or region of interest will not be as informativa as OF .?

There are two versions of optical flow(OF): Feature based (sparse) or dense. In the dense version OF is applied to all the image pixels, while in the sparse one, only certain characteristic feature points are tracked. However, both approaches depend on the tracking of pixel quantities. This is fundamentally different than tracking the whole patch, because in return one obtains a full set of pixel-level correspondences. This is not just a vague statement, but rather carries the idea that the correspondence estimation is not constrained. Remember in the motion case, one would assume that the motion is "rigid" or "articulated" and imposes this prior into the tracking framework. If optical flow is used to estimate the dense trajectories, than this constraint is not assumed and one could as well track the deformable bodies, regardless of the deformation model.

Optical flow provides you more freedom and information about the tracked scene. Yet, my experience is that the motion tracking methods are more robust and reliable. But keep in mind that, by applying human models or temporal models on top of such tracked trajectories (from OF) researchers are capable of developing robust articulated / deformable / model-based tracking algorithms.


Analog Integrator Sampling Rate What would be the effective sampling rate of an op amp analog integrator (such as the one described here ) taking in a rapidly varying signal be? Ie what time resolution would such an integrator have?

The need for this to solve a practical problem has passed, as a different method of testing was ultimately selected, and so I intend this question to be as theoretical as possible, but in order to understand situation better I shall provide the original motivation. I was attempting to use a load cell in order to measure the momentum change of an object colliding with my load cell by taking the force reading of output of the load cell over time and numerically integrating the result. Unfortunately, my data logger only sampled at 1Khz and the impulse event was estimated to take only a few ms, thus providing much too small time resolution to accurately integrate afterwards. So, My thought was to hook the load cell output directly to an analog integrator and I should get a step function output that would read off to me a value proportional to the momentum change. From the transfer function of 1/s, it would seem as though I could do it for any signal regardless of how rapidly changing the signal is, however there must be some practical limit on how fast the integrator can respond. Can I assume a transistor switching speed (on the order of many Mhz and thus much faster than needed here) for the time resolution?

Sampling rate is only defined within discrete-time systems. So, there's really no sampling-rate equivalent that your integrator has.

What your integrator does have is finite bandwidth. It, for example, will not be sensitive for voltage variations at 2.4GHz.

The bandwidth off that circuit involves the characteristics induced by modeling the circuit based on ideal components, but also on limitations such as unintended (parasitic) inductance and capacitance of components and leads, and, most importantly: by the bandwidth of the opamp itself.

You can buy opamps with only a few Hz of bandwidth, you can buy them with multiple GHz of bandwidth, and very many values in between. Each of these have different strengths and disadvantages, so for a measurement system, choose wisely, based on an understanding of the requirements of your measurement system and signal, and at least a solid understanding of electronics and linear electrical networks; I'm not convinced designing a measurement integrator is trivial.

Generally, you might want to take a step back and actually ask over at electronics.SE how to measure a specific signal, describing exactly and in detail what you are observing, rather than assuming an integrator is the solution here. It might be - but it's equally likely you're not on the right track.


FFT of input length 1536 Does anyone knows can i find a FFT of 1536 length input. Its a specification given in 3gpp Lte and we need a transform of 1536 input size which is neither a power of any number i would say. I just need a theoretical idea.

Note that $1536 = 3\cdot 512$. You can use a decimation in time method. See http://www.altera.com/literature/an/an480.pdf

I'll add some extra information here in response to your comment to sansuiso's answer (even though this comment actually referred to the pdf file I mentioned above in this answer). As explained in the above document, you need to compute 3 FFTs of length 512. You do this by taking time samples 0,3,6,9, ... and compute their FFT, then you take the samples 1,4,7,10,... and compute another FFT. Finally you compute the FFT of samples 2,5,8,11,... Then you combine the 3 FFTs using the appropriate twiddle factors to obtain the desired 1536-point FFT. This last step is the radix-3 stage. Have a look at Equation (2) on page 3 of the document (also Fig. 1 on page 4).

Also check out this blog post. It also includes a little Matlab program.